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1.  INTRODUCTION 


Interplanetaty  missions  will  need  rockets  with  high  specific  impulse  and  high  specific  power. 
Nuclear  fusion,  the  power  source  of  the  Sun,  could  power  a  rocket  for  these  interplanetary 
missions.  [l]-[8]  While  a  voyage  to  Mars,  for  example,  would  take  almost  a  year  with 
chemical  propulsion,  a  fusion-powered  rocket  could  get  there  in  months.  The  high 
temperature  of  fusion  plasmas  (already  found  in  eiqperimental  fusion  reactor  concepts)  would 
give  Isp'S  of  lOOO’s  of  seconds.  The  large  energy  release  from  fusion  reactions  would  give 
the  necessary  i^ecific  power. 

If  the  rocket  could  pick  up  more  propellant  when  it  arrived  at  Mars,  it  would  not  need  to 
carry  out  propellant  for  the  voyage  back.  This  would  reduce  the  overall  mass  of  the  inflight 
rocket,  increase  the  thrust-to-weight  ratio,  and  reduce  transit  time.  A  source  of  such 
propellant  is  the  Martian  atmosphere,  which  is  95%  CO2. 

The  JANAF  Thermodynamic  Thbles  [9]  give,  up  to  6000K,  the  enthalpy,  entropy,  and  heat 
capacity  of  several  elements  and  molecules.  However,  fusion  plasmas  have  temperatures  in 
the  lOO’s  of  keV  (10®  K).  By  rederiving  theoretically,  as  the  JANAF  tables  do,  the 
theromodynamic  properties  of  the  fusion  rocket  propellant,  but  using  extended  excitation 
levels,  we  can  estimate  rocket  performance  above  6000  K. 

At  high  temperatures,  as  molecules  in  the  gas  dissociate  into  their  atomic  components, 
chemistry  simplifies.  Although  probability  increases  with  temperature  that  a  molecule  will 
occupy  a  higher  excited  state,  the  number  of  molecules  decreases  with  temperature.  Once 
most  of  the  molecules  have  dissociated,  the  remaining  molecules  have  little  effect  on  the 
thermodynamic  properties  of  the  total  gas. 

At  temperatures  below  6000  K,  the  number  of  chemicals  to  consider  is  overwhelming. 
Computer  codes  exist,  such  as  Selph’s  “ISP  Code”[10],  that  first  identity  chemicals  present 

at  a  given  temperature  and  pressure,  and  then  interpolate  JANAF  table  data  to  find  their 
partial  densities  at  equilibrium.  From  this  they  calculate  the  total  gas  thermodynamic 
properties. 

The  code  described  here  considers:  mixture  of  argon  (Ar),  carbon  (C),  hydrogen  (H), 
nitrogen  (N),  and  oxygen  (O),  up  through  Aiv,  CV,  H//,  N/k,  and  CV;  any  diatomic 
combinations  of  these  elements,  both  neutral  and  singly  ionized.  The  user  can  enter  in 
spectroscopic  data  for  his  own  elements.  Program  results  agree  with  Selph’s  code  at  6000 
K,  for  dominant  ^ecies,  output  graphs  of  nitrogen  species  densities  ( 3000  K  to  30,000  K ) 
and  air  species  mole  fractions  ( 3000  K  to  10,000  K )  match  Cambel[lll. 

The  second  edition  also  has  the  additional  feature  of  being  able  to  calculate  rocket 
performance  based  upon  the  input  of  a  specific  quantity  of  energy  when  running  the  code. 
This  is  e^ecially  valuable  for  calculating  the  high  temperature  performance  of  laser,  fusion, 
and  antiproton  thermal  propulsion  system  which  add  heat  to  a  working  fluid  or  propellant. 


2  THEORY 


2.1  Saha  Equation 

Mass  Action.  Assume  diatomic  molecules,  Aj,  which  dissociate  into  atoms,  A.  At  equi¬ 
librium,  write  the  dissociation  reaction 


A2"2A-  ud,a2, 

(1) 

where  represents  the  dissociation  energy  of  the  reaction. 

According  to  the  law  of  mass 

action  [12],  write  the  equilibrium  constant  of  the  reaction  in  terms  of  the  inventories,  N,  of 

the  individual  species. 

(2) 

Partition  Function.  The  total  energy  of  a  particle  consists  of  its  translational  and 

internal  energies.  Neglect  nuclear  excitation,  and  assume  separable  rotational,  vibrational. 

and  electronic  energies,  to  write  the  particle  Hamiltonian, 

y-tot  —  "Utr  +  "Hrot  +  "Hvib  +  ‘He/. 

(3) 

For  a  particle  partition  function,  write[13] 

Ztot  —  ZtrZfoiZ^ibZtlf 

(4) 

where 

(5) 

(6) 

T 

Zrot  ~  Q.  > 

(7) 

T 

Zvib  «  s — * 

^vib 

(8) 

and  m  represents  the  particle  mass;  h  represents  Planck’s  constant;  k  represents  Boltzmann’s 
constant;  T  represents  the  gas  temperature;  V  represents  the  gas  volume;  «,■  represents 
the  energy  of  the  ith  electronic  excitation  state;  </,  =  2J  +  1  represents  the  degeneracy 
(multiplicity)  of  the  ith  electronic  excitation  state;  J  represents  the  total  angular  momentum 
quantum  number;  0^1  represents  the  characteristic  temperature  for  rotation  of  the  ground 
state  (electronic)  of  the  molecule,  and 


e, 


rol 


Sv^IkT 


(9) 
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where  I  represents  the  moment  of  inertia  of  the  molecule;  represents  the  characteristic 
temperature  for  vibration  of  the  ground  state  (electronic)  of  the  molecule,  and 

=  y,  (10) 

where  u  represents  the  vibration  frequency  of  the  molecule. 

Electronic  excitation  energies,  and  the  degeneracies,  gi,  come  from  spectroscopic 
data[14].  Similarly,  find  Qrot  and  for  diatomic  molecules  in  Huber  k  Hfirzberg[l6]. 
Atomic  species  have  no  vibrational  or  rotational  energy,  |o  thftt 


Zrot  =  Zvih  =  1  (atomic  species). 


ill) 


In  terms  of  the  total  partition  function,  write 


the  reaction  equilibrium  constant  in  Eq.(2) 


Ideal  Gas.  Assume  an  ideal  gas  (a  better  appro.ximation  as  more  molecules  dissociate), 


PV  =  NkT, 


and  substitute  V  into  Eq.(5); 


7  - 

)  P  ' 


where  N  represents  the  total  number  of  particles  in  the  gas,  specifically  “heavy”  particles, 
Nh,  (molecules,  atoms,  ions)  and  electrons, 

+  K.  (15) 

Normalize  all  species  particle  inventories  to  the  total  number  of  nuclei,  A^oi  present  in  the 
gas,  and  define  the  nuclear  fractions  (similar  to  Cambel[ll]): 

VA,  =  (16) 

/  1  7\ 


“  No' 

Nh  , 
y.  =  Yo' 

Ne 

=  lo¬ 


using  Eqs.(lG)-(19)  ill  Eqs.{2)  and  (12)  gives  the  first  Saha  equation: 


VA^iVk  +  Ve) 
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using  the  component  partition  functions  for  the  total  partition  lunctions,  and  simplifying 
the  translational  partition  function  for  and  A.  Note  that  by  taking 
the  second  right-hand-side  term  (in  parentheses)  becomes  the  reduced  mass  of  the  particles. 
Also  note  that,  for  a  given  temperature  and  pressure,  the  value  of  the  right-hana-side  of 
Eq.(20)  is  known  (abbreviated  R). 

Ionization.  Now  consider  ionization  as  well  as  dissociation.  At  equilibrium, 

yl+i +  ,  (21) 

where  j  +  1  represents  the  degree  of  ionization  and  u/j  represents  the  energy  of  the  jth 
ionization.  Equate  the  two  forms  of  the  reaction  equilibrium  constant  to  write 

yA^>{yH  +  y.)  \  h?  )  P  ^  ^ 

Chemical  Reactions.  Now  consider  a  chemical  reaction  between  elements  A  and  B  that 
forms  the  compound  AD.  At  equilibrium, 

AB  ^  A-^  D  —  uii^ab,  (23) 

and  so 

_ yAVB  _  Zel,AZel,B  f-UR,AD\ 

yABiyh+ye)  KmA+rUBJ  P  {ZelZrolZvibjAS  kT  )' 

If  the  compound  AB  ionizes. 


AD  ^  AD+  -h 

e  “  u/./iCt 

(25) 

and 

yAD^Vc 

{kT)^^"^  {ZelZrotZ^ib)AD+  f-Ul,AB\ 

(26) 

VABivk  +  yc) 

V  /i2  ;  p 

{Z,lZrotZ.,b)AB  \  kT  )' 

Saha  Equations.  Solve 

the  following  set  of  non-linear,  coupled,  algebraic  equations: 

Va^Vc 

=  PaIi 

(27) 

yAiyh  +  Vc) 

yA^^Ve 

II 

(28) 

VA+il/A  +  ye) 

VA+^Ve. 

=  Pa2, 

(29) 

yA+Ayh  +  Ve) 

yA+*ye 

= 

(30) 

yA+»(yH-f  y«) 
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BA 

VAtiVk  +  Ve) 

=  ^A8* 

(31) 

=  Rai^ 

(32) 

VA^yh+Vc) 

yB+ye 

—  Rb\, 

(33) 

yaiVk  +  Vc) 

yB^*ye 

=■  RB2f 

(34) 

VB+iyk  +  yt) 

yB**y* 

—  i?B3. 

(35) 

yB+*(yfc  +  ye) 

yB+«y« 

=  Rb4, 

(36) 

VB+^iVk+Ve) 

vl 

=  Rb6, 

(37) 

yB,(y/,  +  y«) 

yBi 

=  i?07, 

(38) 

yBiiyh-^ye) 

yAyu 

=  RAO\y 

(39) 

VAuiyhAye) 

Vad* 

—  RaB2’> 

(40) 

yAoivk  +  ye) 

yh  =  yA  +  va*  +  yA+>  +  yA+»  +  yA+«  +  y^a  +  y^j  + 

(41) 

ys  +  yB+  4"  yB+*  +  yB+>  +  yB+«  +  VBa  +  yfl+  + 

yAB  +  yAB+ ) 

ye  =  yA+  +  2yA+»  +  3y4+a 

+  4yA+«  +  yA+  + 

(42) 

yj5+  +  2i/fl+3  +  3yB+s  +  4yB+4  +  y^  + 
Vab*, 


V^,A  =  yi4  +  yy»+  +  y^+j  +  yA+*  +  yA+«  +  Zy^,  +  2y^+  +  vab  +  Vab^  i  (43) 

yo.B  =  ys  +  yB+  +  ye+*  +  yB+»  +  vb^*  +  Zy^?,  +  2y0+  +  %jad  +  yAB+ ,  (44) 

1  =  yo,A  +  yo.B,  (45) 


where  the  R  (for  “RHS” ,  right-hand-side)  represent  the  known  functions  of  temperature  and 
pressure,  and  yo,A  and  yo.a  represent  the  stoichiometric  fractions  of  elements  A  and  B. 
With  solved  Eqs.(27)-(45),  the  mole  fraction  of  the  kth  species  becomes 


^  =  ^pyt  _  yt 
^  ^oiyh  +  yr)  {yh  +  Vt)' 


Similarly,  the  partial  density  of  the  fcth  species  becomes 


N, 


=  NoVk 


yk 


^Qivh  +  ye)kT  {xjh  +  rjt)  kT ' 


(46) 


(47) 
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2.2  Thermodynamic  Functions 

Sum  the  energies  in  translation,  excitation,  vibration,  electronic  excitation,  dissociation, 
ionization,  and  chemical  reactions  to  get  the  total  energy,  U,  of  the  gas.  Thus, 

U  =  t/(r  +  Urot  +  ^vib  +  ^tl  +  +  Uj  +  U  R.  (48) 


Sum  the  internal  energy  of  the  gas  and  the  thermal  excitation  to  get  the  enthalpy  of  the  gas, 
H: 


^  -U-^PV. 


(49) 


Divide  H  \>y  n  —  NfN^v,  the  number  of  moles,  where  is  Avagadro’s  number,  to  get  the 
molar  enthalpy,  h: 


h 


H 


Av 


N  ’ 


(50) 


•suh.stiting  PV  —  NkT  iuto  the  second  term.  Divide  the  gas  energy  liy  the  number  of  particles 
to  get  the  average  particle  energy. 


u 


£ 

~  N'  _ 

—  Utr  +  +  Urol  +  5e/  +  +  U/  +  tift. 


(51) 


Examine  each  of  the  average  energies  separately. 

Each  of  the  three  degrees  of  freedom  of  translation  has  /cT/2,  and  so  the  average  trans¬ 
lational  energy,  Uir,  becomes 

fi<r  =  \kT.  (52) 

Assume  that  the  average  rotationeil  amd  vibrational  energies  for  diatomic  molecules  have 
kT  each,  a  valid  approximation  at  high  temperatures. 


Urol  ^vib  —  kT . 

Find  the  average  electronic  excitation  energy,  Uei,  from 

-  UigiCxpj-UifkT) 
E.'5.exp(-Ui/A:T) 

_  S.«.gtexp(-u,/fcr) 

Zr, 

Find  the  average  dissociation,  ionization,  and  reaction  energies: 


(53) 


(54) 


— 

=  U/,  and 

=  U/t. 


ud 

U/ 

Ur 
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(55) 

(56) 

(57) 


Recognize  that  Nav^  =  R  —  8.314  J/mol,  the  gas  constant,  so  the  enthalpy  becomes 


A  =  +  +  (58) 

Atomic  species  do  not  have  the  inner  terms  in  parentheses,  since  they  have  neither  rotational 
nor  vibration  energies,  nor  have  reacted  chemically  to  form  a  compound. 

Define  the  constant  pressure  heat  capacity,  cj,, 

(59) 

=  ^  (2  **■ 

Recall  the  definition  of  u,/  from  Eq.{55),  to  find  the  rightmost  term: 


d  _ 
dT^‘‘ 


—  (  exp(-n./A:r)  \ 

ary  z,i  )' 

H  exp{-Ui/kT)  - 

S:,  tiig,exp(-ti,/fcT)  d  ^ 

Zl  dT 

J_ 

kT^  V 

E,-  exp{-Ui/kT)  £,■  Ujgi  exp{-UilkT) 

Zel  Ztl 


Therefore,  write  the  heat  capacity 


Express  a  difference  in  entropy  by[16) 

fr^  dr  P2 

s„-Si=  . 

Operate  on  Cp,  in  Eq.(61),  to  evaluate  the  first  right-hand-side  term: 


d  _  dT 
‘^dr"’’  T ' 


(60) 


(61) 


(62) 


(63) 
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Note  that 


So 


i: 


T*  d  _  dT 
T 


-L 


«»  d 

«,  It 

£fi 

T 


«w» 


a(i/r) 


= 


d 


d{\JT)  ^ exp(-u, /fcT)m, 
=  E^*exp{-a</fcr) 

k 

* 

—^eiZa 


53  «.<7i  exp(-Ui/fcT)m, 


/  (t)  *  /  2*,  5(1/T)  (t)  ’ 

-.-kj 


Substitute  all  the  terms  back  into  the  original  equation  to  get 
S2  -  Si  =  H  |(|  +  (2))  In  ^ 


fiel  ,  ,  Zel,2  ,  Pz 


(64) 


(05) 


(66) 


(67) 


For  convenience,  take  state  1  at  standard  temperature  and  pressure  (STP:  T  =  298.15  K, 
P  =  l  atm). 


Gas  Thermodynamic  Properties.  Equations  (58), (61),  and  (67)  give  formulae  for  the 
enthalpy,  heat  capacity,  and  entropy  difference  for  individual  species  present  in  the  gas.  To 
get  the  entropy,  heat  capacity,  and  entropy  for  the  entire  gas,  sum  the  contributions  from 
each  species.  For  example,  write  the  gas  enthalpy 

where  hi  represents  the  enthalpy  of  the  ith  species,  and  y,7(yfc  4-  ye)  represents  the  mole 
fraction  of  the  ith  species.  Find  gas  heat  capacity  and  enthalpy  difference  in  a  similar 
manner. 


Absorbed  Thermal  Power  To  calculate  the  maximum  absorbed  thermal  power,  assume 
that  the  working  fluid  (propellant)  initially  has  no  enthalpy,  so  Uq  =  0.  Tlic  maximum 
absorbed  thermal  power,  Paba,  becomes 

Paha  =  hlhgaa,  (69) 

where  m  represents  the  mass  flow  rate  of  the  propellant. 
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2.3  Rocket  Applications 

To  estimate  rocket  performance,  make  several  assumptions  [18]: 

1.  Homogeneous  propellant  docs  not  vary  in  composition  throughout  the  rocket  chamber 
and  nozzle. 

2.  The  propellant  obeys  the  perfect  gas  laws. 

3.  The  propellant  flows  adiabatically. 

4.  The  propellant  flows  steadily;  expansion  occurs  isentropically,  without  shocks. 

5.  Exhaust  gases  have  axially  directed  velocity. 

6.  Nozzle  expansion  does  not  change  the  chemical  equilibrium  established  in  the  chamber 
(“frozen  flow”). 

7.  Propellant  contains  gaseous  species  only. 

The  diatomic  species  do  not  exacly  follow  the  ideal  gas  assumptions,  since  they  have 
degrees  of  freedom  in  rotation  and  vibration.  At  high  temperatures,  as  the  diatomic  species 
dissociate,  the  ideal  gas  assumption  improves.  As  the  nozzle  expands  and  the  temperature 
and  pressure  decrease,  the  dissociation/ionization  equilibrium  changes.  However,  assume 
that  expansion  occurs  quickly  enough  that  the  equilibrium  at  the  exit  remains  unchanged 
from  the  chamber,  and  that  the  expansion  occurs  isentropically  (reversibly).  Use  the  frozen 
flow  approximation,  rather  than  calculate  the  equilibrium  conditions  at  both  the  chamber 
and  exit,  since  the  Saha  equations  do  not  give  correct  answers  at  very  low  pressures.' 

Since  the  expansion  occurs  isentropically, 


where  7  =  5/3  =  <^/c„  for  an  ideal  gas.  The  heat  capacity  for  an  ideal  gas  remains  constant 
at  Cp  =  /?5/2,  so  conservation  of  energy  yields 

^(<4-"’)  =  K.-h.  (71) 

=  c,{T„-T). 

Assume  that  the  axial  velocity  of  the  propellant  in  the  chamber,  v,  is  small  compared  to 
the  exit  velocity,  v*,,  to  solve  for  the  exit  velocity; 

=  y2Cp(Te,-“r,  (72) 

‘The  Saha  equations  describe  equilibrium  brought  about  from  collisions.  At  low  pressures,  fewer  collisions 
occur,  and  cquilbrium  occurs  by  way  of  iongcr-range  interactions  (“coronal  equilibrium”).  This  project 
considered  only  Saha  equilibrium. 


Convert  R  from  J/K/mol  to  J/K/kg  by  dividing  R  by  the  average  molecular  weight,  JTW, 

(Vh  +  ye) 


The  exit  velocity  gives  thrust,  F,  specific  impulse,  I^p,  and  jet  power,  Pjet'. 


(74) 


F 

I»p 

Piet 


mva, 
t'e.  . 


,  and 


(75) 

(76) 

(77) 
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3  USER’S  MANUAL 


This  section  tells  how  to  use  the  code,  what  the  input  parameters  mean,  what  techniques  work 
best  for  the  desired  results,  where  to  look  for  output,  and  how  to  create  new  input  elements.  If 
the  user  has  trouble,  the  author  receives  internet  email  at  the  address  nachtriebCuiuc .  edu, 
or  US  mail  care  of  the  Phillips  Laboratory. 

3.1  Starting  Up 

By  making  a  directory  on  the  hard-drive  to  which  all  the  program  files  may  be  copied,  the 
user  can  keep  input  and  output  files  separate  from  other  programs,  and  can  keep  track  of 
the  output  files  the  program  creates.  To  do  so,  at  the  prompt  (something  like  C:\>),  type; 

C:\>  md  saha 

C:\>  copy  a;*.*  c:\saha 

C:\>  cd  \saha 

C:\SAHA> 

These  commands  make  a  directory  named  SAHA,  copy  the  program  files  from  floppy  to  the 
new  directory,  and  then  change  into  the  new  directory. 

To  run  the  code,  type  the  program  name,  ab: 

C;\SAHA>  ab 

Several  examples  might  prove  helpful  for  demonstrating  how  to  operate  the  code. 

3.2  Example:  Carbon  Dioxide,  One  Temperature 

A  few  lines  of  introduction  greet  the  user,  after  which  the  program  prompts  the  user  to 
indicate  which  two  elements  he  wishes  the  program  to  work  with. 

Program;  AB.FOR,  ver.  1.3;  Author;  R.Nachtrieb,  Jan  93; 

OLAC-PL/RKFE  Edwards  AFB  CA  93523-5000 
internet  email:  nachtriebfluiuc.edu 


Ar,  C,  H,  N,  0 

Which  two  (2)  elements  to  use? 

[C  .0  ] 

The  user  types  the  letters  corresponding  to  the  symbols  of  the  elements  he  wishes  to 
consider.  The  program  suggests  four  elements:  carbon,  hydrogen,  nitrogen,  and  oxygen, 
as  indicated  by  C,  H,  N,  0.  In  hard  brackets,  the  program  displays  the  current  (default) 
setting.  If  the  user  wants  to  use  the  default  settings,  he  simply  presses  <Enter>.  Otherwise, 


he  can  type  in  two  element  symbols,  separated  by  either  a  space  or  a  comma.  Assume  that 
the  user  accepts  the  default  elements  and  presses  <Enter>. 

Once  the  user  has  entered  in  his  selections,  the  program  echoes  back  the  selection.  Note 
that  the  program  is  case  sensitive,  so,  for  example,  it  does  not  recognize  n  as  nitrogen;  the 
user  should  use  N  for  nitrogen. 

Using  C  ,0  ,C0 

Stoichiometric  Proportions  of  C  ,0  ? 

[  l.OOOE+00,  2.000E+00] 

The  program  then  asks  the  user  for  the  relative  proportions  of  each  element.  For  example, 
if  the  user  wished  to  analyze  carbon  dioxide,  the  user  could  indicate  that  there  was  one  part 
of  carbon  for  two  parts  of  oxygen.  Note  that  this  program  docs  not  consider  polyatomic 
molecules,  only  diatomic;  but  the  user  can  approximate  gases  like  CO2  by  giving  the  same 
stoichiometric  proportions  as  CO2. 

Proportions  of  C  ,0  l.OOOE+00  2.000E+00 

Chamber,  exit  pressures  (atm)? 

C  i.oooE+00.  o.oooE+003 

The  program  prompts  the  user  for  the  chamber  pressure,  in  atmospheres,  to  use  in  the 
Saha  equations.  The  exit  pressure  is  used  to  calculate  ideal  rocket  performance.  At  low 
pressures,  the  Saha  equations  no  longer  give  the  correct  equilibrium  species  distribution, 
since  collisions  between  gas  molecules  occur  infrequently.  Thus,  while  the  program  may  give 
answers  for  chamber  pressures  below  10'^  atm,  the  user  should  not  trust  them.  The  program 
does  not  invoke  the  Saha  equations  to  calculate  equilibrium  conditions  at  the  nozzle  exit,  so 
the  user  may  enter  zero  pressure  (vacuum)  for  the  exit. 

pressl,press2  (atm)  l.OOOE+00  O.OOOE+00 

Specify  Temperature  (T)  or  Absorbed  Power  (P)? 

[T3 


The  program  asks  the  user  either  to  specify  the  chamber  temperature  and  calculate  the 
absorbed  power,  or  to  specify  the  absorbed  power  and  calculate  the  chamber  temperature. 
The  program  actually  uses  the  chamber  temperature  as  an  independent  variable  and  calcu¬ 
lates  the  absorbed  power.  When  the  users  specifies  a  desired  absorbed  power,  the  program 
varies  the  chamber  temperature  until  the  calculated  absorbed  power  matches  the  desired  ab¬ 
sorbed  power.  A  later  example  describes  specifying  the  absorbed  power.  For  this  example, 
choose  to  specify  the  chamber  temperature,  T. 
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T 

Chamber  temperature  range:  tempi (K), temp2CK) .howmany? 

C  6.000E+03,  3.000E+04,  1] 

The  prograjn  asks  the  user  for  the  temperature  range  to  consider.  If  the  user  wishes  to 
consider  a  single  temperature,  he  types  in  the  low  and  high  temperatures  of  the  temperature 
range,  and  then  types  a  1.  The  program  checks  to  see  how  many  temperatures  the  user 
wants  to  consider:  if  the  user  has  entered  in  a  1,  the  program  ignores  the  second  (high) 
temperature  and  considers  only  the  first.  Thus,  to  consider  a  gas  with  temperature  6000  K, 
the  user  types  6000,  followed  by  a  space,  another  temperature  (any  number  will  do,  since 
the  program  ignores  it),  and  then  a  1. 

The  program  can  run  for  temperatures  below  6000  K.  However,  since  it  does  not  consider 
polyatomic  molecules,  the  program  looses  accuracy  the  further  below  6000  K  the  user  specifies 
the  temperature.  Below  6000  K,  the  user  should  use  Selph’s  “Isp”  code. 

tempi (K)  .temp2(K)  .howmany  6.000E-'-03  3.000E+04  1 

View  loop  progress:  inner, middle, outer? 

[  20.  20,  20] 

The  numbers  that  the  user  enters  tell  the  program  how  many  loop  iterations  the  program 
should  count  before  printing  a  progress  message  to  the  screen.  On  personal  computers,  the 
Saha  equations  may  take  a  long  time  to  converge,  especially  if  the  computer  does  not  have 
a  numeric  co-processor.  To  reassure  the  user  that  the  program  runs  properly,  the  program 
prints  convergance  factors  to  the  screen,  with  the  number  of  times  the  program  loops  has 
iterated.  Using  fast  computers,  printing  out  to  the  screen  every  loop  slows  down  the  program; 
but  using  slow  computers,  the  user  can’t  tell  if  the  program  works  or  not,  since  the  program 
taJees  a  long  time  to  show  signs  of  activity. 

Sugestions.  The  firet  time  the  user  runs  the  program,  try  entering  20  20  20.  If  the 
program  seems  to  iterate  the  numbers  as  fast  as  it  can,  that  is,  the  program  does  not  pause 
between  progress  updates,  try  increasing  the  numbers  to,  say,  100  100  20.  The  first  two 
numbers  need  not  be  greater  than  100,  since  the  loops  usually  converge  in  fewer  than  100 
intcrations. 

The  program  is  Mono”  when  the  final  residual  (sum  of  all  the  absolute  v-.iluos  of  the 
differences  between  the  left-hand-sides  and  the  right-hand-sides  of  the  Saha  equations)  drops 
below  c,  where  e  =  10"^.  When  the  outer  loop  residual  gets  down  around  10“®,  it  starts 
bouncing  around.  If  the  residual  doesn’t  drop  below  e  within  a  specified  number  of  outer 
loop  iterations  (the  last  of  the  three  numbers  the  user  types),  the  program  stops.  Usually 
20  outer  loop  iterations  suffices. 

Note  that  only  the  final  of  the  three  numbers  affects  the  accuracy  of  the  program  results. 
In  all  cases,  the  loop  residuals  must  drop  below  e.  However,  if  the  program  can’t  quite  get 
the  outer  loop  residual  below  e,  but  comes  close,  the  user  can  stop  calculations  after,  say,  20 
iterations. 
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seein.seemid.maxout  20  20  20 


mass  flow  rate  (kg/s)? 
£  l.OOOE+00] 


Finally,  the  program  asks  the  user  for  the  mass  flow  rate.  The  program  uses  the  mass 
flow  rate  to  determine  the  thrust  an  ideal  rocket  nozzle  produces. 

After  reading  the  mass  flow  rate,  the  program  has  all  the  information  it  needs  to  run  a 
single  temperature  case.  It  echoes  back  the  information,  and  then  starts  to  work.  First  it 
reads  in  spectroscopic  data  from  disk. 

mdot :  1 


Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 
Getting  data  from 


file  C.O.in 
file  C.l.in 
file  C_2.in 
file  C_3.in 
file  C_4.in 
file  0_0.in 
file  O.l.in 
file  0_2.in 
file  0_3.in 
file  0_4.in 


Then  the  program  starts  to  work  on  solving  the  Saha  equations.  It  writes  to  the  screen 
what  temperature  (in  K)  and  pressure  (now  in  pascals,  P;  it  considers.  It  also  writes  six 
columns  of  numbers,  which  tell  the  number  of  loop  iterations  and  the  residuals  for  the  three 
program  loops. 


temp(K),  press 1 (Pa)  6.000E+03  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyk-l  out. loop  residual 

2.000E+01  2.546E-07 


When  the  program  has  linished  with  a  temperature,  it  prints  out  the  maximum  absorbed 
thermal  power  and  the  jet  power,  thrust,  and  specific  impulse  that  an  ideal  rocket  would 
produce.  The  program  also  tells  which  file  to  examine  to  get  the  complete  information  it 
calculated. 


Pabs,  Pjet  (W)  2.541E+06  5.556E+06 

Thrust (N),Iap(s)  3.333E+03  3.399E+02 

Look  for  output  in  file  (CO. out  ) 

Press  <Enter>  to  exit  program. . . 
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Note  that  for  the  last  example,  the  jet  power  exceeds  the  maximum  absorbed  thermal 
power.  This  counterintuitive  result  comes  from  assumptions  made  for  the  ideal  rocket.  At 
higher  chamber  temperatures  the  maximum  absorbed  thermal  power  exceeds  the  jet  power, 
as  it  should. 

Single- Temperature  Output  File.  We  can  print  a  hardcopy  of  the  output  file,  CO .  out, 

by  typing  at  the  MS-DOS  prompt: 

C:\SAHA>  type  co.out  >  prn 

However,  after  several  runs,  this  may  amount  to  a  large  stack  of  paper.  An  easy  way  to  view 
the  program  on  the  screen  uses  the  “Browse”  function  of  Frailey’s  DCOM  program  (or  an 
ASCII  text  editor). 

The  first  few  lines  of  the  output  file  echo  back  the  users  input. 

Using  C  ,0  ,CQ 

Stoichiometric  Proportions  of  C  ,0 
[  3.333E-01.  6.667E-01] 

chamber,  exit  pressures  (atm) 

C  1.013E+05,  O.OOOE+00] 

Temperature  range:  templ(K) ,temp2(K) .howmany 
[  6.000E+03,  3.000E+04.  1] 

mass  flow  rate  (kg/s) 

C  l.OOOE+00] 

Next,  the  file  contains  the  partition  functions  for  all  species  considered  in  the  program. 
Partition  functions  for  atomic  species  appear  in  the  first  five  columns,  in  increasing  order 
of  ionization.  The  sixth  and  seventh  columns  contain  partition  functions  for  homonuclear 
diatomic  species,  also  in  order  of  increasing  charge.  The  three  rows  of  data  contain  the 
electronic  (labeled  el  at  the  far  left),  rotational  (rot),  and  vibrational  (vib)  partition  func¬ 
tions,  respectively.  (Note  that  the  rotational  and  vibrational  partition  functions  for  the 
atomic  species  have  1  .OOOE+00.) 

Partition  function  data  for  the  heteronuclear,  diatomic  molecule  follow  the  elemental 
l)aitition  function  data.  Only  the  diatomic  columns  contain  data,  in  the  same  electronic, 
rotational,  and  vibrational  row  order. 


Partition  Function:  C 


atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic, 

diatomic, 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

9.377E+00 

5.939E+00 

l.OOOE+00 

2.000E+00 

l.OOOE+00 

8.794E+01 

5.570E+01 

vib 

l.OOOE+00 

l.OOOE+00 

1 . OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.292E+03 

2.515E+03 

rot 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.249E+00 

3.090E+00 
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Partition  Function:  0 


atomic , 

atomic . 

atomic. 

atomic. 

atomic, 

diatomic, 

diatomic, 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el  8.947E+00 

4.016E+00 

8.604E+00 

5.646E+00 

l.OOOE+00 

8.005E+01 

3.594E+01 

vib  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.886E+03 

2.466E+03 

rot  l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

2.640E+00 

2.190E+00 

P^o:tition  Function:  CO 

atomic. 

atomic. 

atomic. 

atomic , 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

8.390E+01 

4.540E+01 

vib 

3.227E+01 

1 . 073E+02 

rot 

2.844E+00 

3.185E+03 

The  output  file  contains  the  right-hand-sides  of  the  Saha  equations,  found  in  Eqs.(27)- 
(45). 

Right-hand-sides  of  Saha  Equation 


atomic. 

atomic. 

— 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C  1.998E-07 

S.122E-19 

O.OOOE+00 

O.OOOE+00 

0 .  OOOE+00 

1 . 240E+00 

6.017E-08 

0  1.502E-09 

5.827E-27 

0 .  OOOE+00 

O.OOOE+00 

O.OOOE+00 

1.083E+01 

2.434E-08 

CO 

1.782E-04 

8.278E-10 

The  solution  to  the  Saha  equations  gives 

the  nuclear  fractions,  that  is,  the  species  inven- 

tory  divided  by  the  total  number  of  nuclei. 

Nuclear  fraction8,y 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C  1.265E-04 

6 . 198E-07 

7.789E-21 

O.OOOE+00 

O.OOOE+00 

9.868E-09 

2.914E-11 

0  3.068E-01 

1.131E-05 

1.616E-27 

O.OOOE+00 

O.OOOE+00 

1.331E-02 

7.944E-06 

CO 

3.332E-01 

6 . 767E-06 

yh,ye,yh+ye,beta  6.535E- 

■01  2.664E-05  6.535E- 

•01  4.076E- 

■05 

Dividing  the  nuclear  fractions  by  the  sum  of  the  electron  fraction  and  the  heavy 

fraction, 

yh+ye,  gives  the  mole  fractions. 

Hole  fractions 

,y/ (yh+ye) 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

C  1.935E-04 

9.484E-07 

1.192E-20 

O.OOOE+00 

O.OOOE+00 

1.510E-08 

4.458E-11 

0  4.695E-01 

1.730E-05 

2.473E-27 

O.OOOE+00 

O.OOOE+00 

2.036E-02 

1.216E-05 

CO 

5.099E-01 

1.036E-05 

e-  4.076E-05 
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• 

The  output  file  contains  the  densities  of  the  species,  per  cubic  meter. 

Species  densities  (a^-S) 

• 

atomic.  atomic. 

atomic. 

atomic. 

atomic , 

diatomic. 

diatomic, 

neutral  +1 

+2 

+3 

+4 

neutral 

+1 

C 

2.367E+20  1.16CE+18 

1.458E+94 

O.OOOE+00 

O.OOOE+00 

1.847E+16 

5.453E+13 

- 

0 

5.743E+23  2.117E+19 

3.025E-03 

O.OOOE+00 

O.OOOE+00 

2.490E+22 

1.487E+19 

CO 

6.236E+23 

1.267E+19 

e- 

4.986E+19 

One  finds  the  molar  heats  of  formation  of  the  species  by  adding  half  the  dissociation 

energy  to  the  ionization  energies,  less  the  reaction  energy.  The  output  file  lists  heats  of 
reaction  for  the  species;  subsequently,  the  enthalpies  listed  include  the  heats  of  reaction  to 

give  “absolute”  enthalpy.  The  output  file  also  lists  the  total  gas  enthalpy. 

Heat  of  formation  (J/mol) 

atomic,  atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral  +1 

+2 

+3 

+4 

neutral 

+1 

C 

2.996E+05  1.299E+06 

3.462E<«-06 

7.710E+06 

1.343E+07 

O.OOOE+00 

i.l67E+06 

0 

2.465E+05  1.455E+06 

4.S74E+06 

9.449E’i’06 

1.632E+07 

O.OOOE+00 

1.158E+06 

CO 

“5.152E+05 

8.356E+05 

• 

Enthalpies  (J/mol) 

atomic,  atomic. 

atomic , 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral  +1 

+2 

+3 

+4 

neutral 

+1 

C 

4.305E+0S  1.512E+06 

3,864E+06 

8.483E-^06 

1.471E+07 

2.369E+05 

1.399E+06 

0 

3.745E+05  1.686E+06 

5.080E+06 

1.038E+07 

1.785E+07 

2.310E+05 

1.387E+06 

CO 

-2.812E+05 

-2.850E+05 

e- 

1.247E+05 

CO 

gas  3.727E+04 

The  output  file  also  contains  species  and  gas  heat  capacities  and  entropy  differences  (from 

STP). 

Heat  capacities  (C_p)  (J/K/mol) 

atomic,  atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral  +1 

+2 

+3 

+4 

neutral 

+1 

C 

2.318E+01  2.084E+01 

2.083E+01 

2.079E+01 

2.079E+01 

4.169E+01 

3.971E+01 

0 

2.227E+01  2.223E+01 

2.166E+01 

2.080E+01 

2.079E+01 

4.024E+01 

4.028E+01 

CO 

4.100E+01 

4.000E+01 

e- 

2.079E+01 

• 

CO 

gas  3.219E+01 

• 
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Entropy  diff  (from  STP)  (J/K/mol) 


atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

C 

2.318E+01 

2.084E+01 

2.083E+01 

2.079E+01 

2.079E+01 

4.169E+01 

0 

2.227E+01 

2.223E+01 

2.166E+01 

2.080E+01 

2.079E+01 

4.024E+01 

CO  4.100E+01 

e-  6.240E+01 
CO  gas  9.313E+01 


diatomic , 
+i 

3.971E+01 

4.028E+01 

4.000E+01 


Finally,  the  output  file  lists  the  rocket  performance. 


Thrust(N),Isp(s)  3.333E+03  3.399E+02 

Pabs,  Pjet  (W)  2.541E+06  5.556E+06 

The  program  stores  the  user’s  input  parameters  in  file  AB. in,  and  offers  the  parameters 
as  defaults  (in  hard  brackets)  the  next  time  the  user  runs  the  program. 


3.3  Example:  Carbon  Dioxide,  Multiple  Temperatures 

To  examine  multiple  temperatures,  the  user  enters  in  data  in  much  the  same  way  as  for  a 
single  temperature.  Once  again,  run  the  program, 

C:\SAHA>  ab 

and  provide  the  program  with  the  information  it  requests: 

Program:  AB.FOR,  ver.  1.3;  Author:  R.Nachtrieb,  Jan  93; 

OLAC-PL/RKFE  Edwards  AFB  CA  93523-5000 
internet  email:  nachtrieb9uiuc.edu 


Ar,  C.  H.  N,  0 

Which  two  (2)  elements  to  use? 

[C  .0  ] 

Using  C  ,0  ,C0 

Stoichiometric  Proportions  of  C  ,0  ? 

[  l.OOOE+00,  2.000E+00] 

Proportions  of  C  ,0  l.OOOE+00  2.000E+00 

Chamber,  exit  pressures  (atm)? 
t  l.OOOE+00,  O.OOOE+00] 
pressl,press2  (atm)  l.OOOE+OO  0,OOOE+00 
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Specify  Temperature  (T)  or  Absorbed  Power  (P)? 

[T] 

T 

Chamber  temperature  range:  templ(K) ,temp2(K) .howmany? 

C  6.000E+03,  3.000E+04,  1] 

Now  the  user,  instead  of  selecting  only  one  temperature,  inputs  the  temperature  range 
he  wishes  to  consider,  and  the  number  of  divisions  within  that  temperature  range.  If  the 
user  enters  in  2  for  howmany  divisions,  the  program  considers  the  lower  and  upper  bounds; 
if  the  user  enters  3,  the  program  considers  the  lower  and  upper  bounds,  and  a  temperature 
midway  between  the  two,  and  so  on. 

Try  entering  in 

6e3  30e3  3 

The  program  echoes  back  the  user’s  selection,  and  then  prompts  for  a  new  input: 
tempi (K) ,temp2(K) .howmany  6.000E+03  3.000E+04  3 

Print  density  (1)  or  mole  fraction  (2)  ? 

[13 

When  the  user  selects  more  than  one  temperature  for  the  program  to  consider,  the  program 
outputs  information  in  a  different  way:  instead  of  writing  all  information  to  one  file,  the 
program  write  specific  information  to  many  files.  This  makes  it  easy  for  the  user  to  import 
a  file  of,  say,  species  densities  vs.  temperature  into  a  graphing  program.  In  this  way,  the 
user  simply  looks  at  the  particular  file  to  see  enthalpy,  heat  capacity,  entropy,  mole  fraction, 
density,  or  rocket  performance.  Unfortunately,  MS-DOS  limits  the  number  of  files  a  user 
(actually,  the  program)  can  have  open  at  one  time.  Thus  the  user  must  choose  whether  to 
make  files  that  contain  species  densities,  or  to  make  files  that  contain  species  mole  fractions. 
Note  that  if  the  user  wants  buth  density  and  mole  fraction,  he  can  run  the  program  twice, 
and  specify  density  the  first  time  and  mole  fraction  the  second  time. 

The  rest  of  the  inputs  follow  the  single  temperature  example. 

de’smol  [1] 


View  loop  progress:  inner .middle, outer? 

[  100,  100,  20] 

seein.seemid.maxout  100  100  20 

mass  flow  rate  (kg/s)? 

[  l.OOOE+00] 
mdot :  1 
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After  the  user  has  input  the  information,  the  program  reads  atomic  spectra  data  from 
input  files,  and  starts  to  work. 

Getting  data  from  file  C_0.in 
Getting  data  from  file  C_l.in 
Getting  data  from  file  C_2.in 
Getting  data  from  file  C_3.in 
Getting  data  from  file  C.4.in 
Getting  data  from  file  0.0. in 
Getting  data  from  file  O.l.in 
Getting  data  from  file  0_2.in 
Getting  data  from  file  0_3.in 
Getting  data  from  file  0_4.in 

Since  the  program  now  considers  multiple  temperatures,  it  gives,  for  each  temperature, 
the  pressure,  progress  updates  for  the  loop  calculations,  and  rocket  performance. 

temp(K),  pressl(Pa)  6.000E+03  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 

2.000E+01  2.546E-07 

Pabs,  Pjet  (W)  2.541E+06  5.556E+06 

Thrust (N),Isp(s)  3.333E+03  3.399E+02 

temp(K),  pressl(Pa)  1.800E+04  1.013E+05 

in. loop  wyhX-yOX  mid. loop  wyh-1  ou  .loop  residual 

2.000E+01  1.126E-06 

Pabs,  Pjet  (W)  7.476E+07  4.561E+07 

Thrust  (N).Isp(s)  9.551E+03  9.740E-r02 

temp(K),  pressl(Pa)  3.000E+04  1.013E+05 

in, loop  wyhX-yOX  mid. loop  wyh-1  out. loop  residual 

2.000E+01  1.859E-06 

Pabs,  Pjet  (W)  i.l66E+08  9.779E+07 

Thrust (N).Isp(s)  1.398E+04  1.426E+03 

Finally,  the  program  tells  the  user  which  files  to  examine  for  the  output  data.  Note 
that  the  filenames  follow  mnemonic  convention:  to  see  a  file  containing  maximum  absorbed 
thermal  power,  jet  power,  rocket  thrust,  and  specific  impulse  as  a  function  of  chamber 
temperature,  see  file  CO.rkt.out;  for  carbon  species  densities,  see  file  C.den.out,  and  so 
on. 

Look  for  output  in  files; 

(C_eth.out  ) 

(O.eth.out  ) 
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(CO.eth.out  ) 

(C.Cp . out  ) 

(O.Cp . out  ) 

(C0_Cp . out  ) 

(C.etr.out  ) 

(0_etr.out  ) 

(CO.etr.out  ) 

(C.den.out  ) 

(O.den.out  ) 

(CO.den.out  ) 

(C0_rkt.out  ) 

Press  <Enter>  to  exit  program. . . 

Multiple-Temperature  Output  Files.  When  the  users  selects  multiple  temperatures, 
the  output  files  created  contain  data  in  columns.  For  example,  the  file  CO_rkt .  out  contains 
rocket  performance  data,  as  a  function  of  temperature: 

CO  rocket  performance  parameters 

temp(K)  Pabs(W)  Pjet(W)  thrust(N)  IspCs'' 

6.000E+03  2.541E+06  5.5S6E+06  3.333E+03  3.399E+02 

1.8C0E+04  7.476E+07  4.56iE+07  9.551E+03  9.740E+02 

3.000E+04  1.166E+08  9.779E+07  1.398E+04  1.426E+03 

The  program  creates  three  output  files  for  density  (or  mole  fraction),  enthalpy,  heat 
capacity,  and  entropy.  Thus,  file  C..den.out  contains  the  desities  of  .all  the  carbon  atomic 
species,  the  homonuclear  diatomic  carbon  species,  and  the  total  electron  density.  Similarly, 
file  0_den.out  contains  density  for  oxygen  species,  and  the  total  electron  density.  The  file 
C0_den .  out  contains  density  of  carbon  monoxide  species,  and  the  total  electron  density.  The 
program  writes  the  same  total  electron  densities  to  each  file  as  a  safety  check:  if  the  user 
plots  all  the  densities  on  the  same  graph,  discrepancies  in  electron  density  alert  the  user  to 
an  error  (for  example,  plotting  a  file  from  a  different  program  run). 

3.4  Example:  Nitrogen,  Multiple  Temperatures 

Suppose  the  user  wants  to  consider  a  single  element,  say  nitrogen.  Run  the  program,  and 
Eifter  it  asks  for  the  elements  to  consider, 

Ar.  C,  H,  N,  0 

Which  two  (2)  elements  to  use? 

[C  .0  ] 

enter  in 

N  N 
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The  program  echoes  back  the  users  selection,  and  asks  for  the  stoichiometric  proportions. 
Using  N  .N  ,m 

StoichioBketric  Proportions  of  N  ,H  ? 

C  l.OOOE+00,  2.000E+00] 

At  this  point,  the  user  can  enter  any  proportions  he  desires,  since  both  species  are  the  same. 
However,  a  much  better  way  exists: 

Hint.  When  considering  only  one  element,  type  0  for  one  of  the  proportions.  This  speeds 
up  the  code  considerably,  since  it  no  longer  has  to  consider  “heteronuclear”  diatomic  species. 
Try 

1  0 

after  which  the  program  responds, 

Proportions  of  N  ,N  O.OOOE+00  l.OOOE+00 
After  this,  the  program  functions  as  before. 

3.5  Example:  Specified  Absorbed  Power 

Now  run  the  program  again,  but  this  time  specify  a  maximum  absorbed  power  instead  of  a 
chamber  temperature. 

Ar,  C,  H,  N.  0 

Which  two  (2)  elements  to  use? 

[C  .0  3 

Using  C  ,0  ,C0 

Stoichiometric  Proportions  of  C  ,0  ? 

C  l.OOOE+00,  2.000E+00] 

Proportions  of  C  ,0  l.OOOE+00  2.000E+00 

Chamber,  exit  pressures  (atm)? 

[  l.OOOE+00,  O.OOOE+00] 
pressl,press2  (atm)  l.OOOE+00  O.OOOE+00 

Specify  Temperature  (T)  or  Absorbed  Power  (P)? 

CT]  P 
P 


The  program  then  does  not  ask  which  temperature  range  to  examine,  but  goes  directly 
to  the  loop  progress  and  mass  flow  rate  inputs. 
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View  loop  progress;  inner, middle, outer? 

[  100,  100,  20] 

seein,8eemid,maxout  100  100  20 


mass  flow  rate  (kg/s)? 
C  l.OOOE+00] 
mdot :  1 


The  program  then  gets  the  input  data  files  and  calculates  the  lower  and  upper  bounds 
on  the  maximum  absorbed  thermal  power. 


Getting 

Getting 

Getting 

Getting 

Getting 

Getting 

Getting 

Getting 

Getting 

Getting 

Getting 


data  from  file  C.O.in 
data  from  file  C.l.in 
data  from  file  C.2.in 
data  from  file  C.3.in 
data  from  file  C_4.in 
data  from  file  0.0. in 
data  from  file  0.1. in 
data  from  file  0.2. in 
data  from  file  0.3. in 
data  from  file  0.4. in 
lower  absorbed-power  bound... 


Getting  upper  absorbed-power  bound... 


Pmin.Pmax  (W)  2.541E+06  1.166E+08 


2.000E+01  2.S46E-07 
2.000E+01  1.859E-06 


The  program  results  have  been  verified  for  pressure  P  =  1  atin  between  temperature 
T  =  6000  K  and  T  =  30000  K.  The  program  thus  takes  T  —  6000  K  and  T  =  30000  K  as 
the  temperature  bounds  and  assumes  valid  results  within  those  bounds.  Using  the  chemical 
composition  of  the  propellant  and  the  mass  flow  rate,  the  program  calculates  the  maximmii 
absorbed  thermal  powers  at  T  =  6000  K  and  T  =  30000  K  and  displays  them  to  the  user. 
The  program  then  asks  the  user  for  the  desired  absorbed  thermal  power,  Pwant.  The  desired 
thermal  power  must  be  between  the  absorbed  thermal  powers  calculated  at  the  boundry 
temperatures. 


What  power  to  absorb?  (W)  le? 


resP 

6.476E+00 
3 . 243E+00 
9.7eiE-01 
2.619E-01 
2.144E-01 


epsP 

l.OOOE-02 

l.OOOE-02 

l.OOOE-02 

l.OOOE-02 

l.OOOE-02 


temp 

1.800E+04 
1.200E+04 
9 . OOOE+03 
7.500E+03 
8.250E+03 


out . loop 

2.000E+01 

2.000E+01 

2.000E+01 

2.000E+01 

2.000E+01 


residual 

1.126E-06 

6.323E-07 

2.923E-07 

3.104E-O7 

2.725E-07 
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6.200E-02  l.OOOE-02  7.875E+03 
6.580E-02  l.OOOE-02  8.063E+03 
6.543E-04  l.OOOE-02  7.969E+03 


2.000E+01 

2.000E+01 

2.000E+01 


1.333E-07 
3.229E-07 
1 . 179E-07 


Psant.Pabs.temp  l.OOOE+07  9.993E+06  7.969E+03 
Pabs,  Pjet  (W)  9.993E+06  7.808E-*-06 

Thrust (N),Isp(8)  3.952E+03  4.030E+02 

The  program  then  searches  for  the  temperature  at  which  the  calculated  absorbed  thermal 
power  matches  the  desired  absorbed  thermal  power,  to  within  1%.  The  program  displays 
progress  of  decrease  in  the  residuals  and  the  search  temperatures.  If  the  user  enters  a  desired 
absorbed  power,  Pwant,  outside  the  boundries  specified  by  the  absorbed  thermal  powers  at 
T  =  6000  K  and  T  =  30000  K,  the  program  returns  the  user  to  the  beginning  of  the  data 
input. 

Since  the  program  iterates  several  temperatures  to  match  the  calculated  absorbed  thermal 
power  to  the  desired  absorbed  thermal  power,  it  asks  for  only  one  desired  power.  However, 
the  absorbed  thermal  power  and  the  chamber  temperature  have  a  one-to-one  correspondence. 
The  user  can  plot  the  absorbed  thermal  power  and  jet  power  as  a  function  of  temperature, 
as  in  Fig.  1  and  find  the  temperature  corresponding  to  the  desired  absorbed  thermal  power. 


Figure  1:  Maximum  Absorbed  Thermal  Power  and  Ideal  Rocket  Jet  Power  vs.  Chamber 
Teniperaturc,  CO2,  P  =  1  atm.  Upper  curve  represents  Absorbed  Thermal  Power. 
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3.6  Example:  New  Element 

Suppose  the  user  wants  to  run  the  program  with  an  element  the  program  doesn’t  yet  recog¬ 
nize.  The  program  provides  a  way  for  the  user  to  enter  the  information  the  program  needs 
to  run,  and  save  that  information  to  the  input  files. 

The  user  runs  the  program,  and  enters  in  the  elements  he  wishes  to  consider,  C  and  Pn.* 

Ar,  C,  H.  N,  0 

Which  tuo  (2)  elements  to  use? 

[H  .N  ]  C  Pn 

Oops!  Couldn’t  find  Pn  in  file  (ABX.in) 

Make  a  new  element?  (l^yes) 

Cl] 

The  file  ABX.in  contains  input  information  about  the  elements  the  program  recognizes,  and 
the  homonuclear  diatomic  molecules  formed  by  those  elements;  specifically,  the  program 
looks  in  ABX .  in  for  the  highest  ionization  state  and  the  mass  of  the  atomic  species,  for  dis¬ 
sociation  and  ionization  energy  of  the  homonuclear  diatomic  molecule,  and  for  characteristic 
temperatures  of  rotation  and  vibration.  The  program  looks  through  the  input  file  ABX.in 
for  the  element  Pn,  but  doesn’t  find  it.  It  then  asks  the  user  if  he  wants  to  define  a  new 
element,  and  provide  its  characteristic  information.  If  the  user  says  no  (0),  the  program  asks 
the  user  for  two  new  elements.  If  the  user  says  yes  (1),  the  program  continues. 

Let’s  make  a  new  element... 

Highest  atomic  state  (l+highest  ionization  state)? 

CO] 

The  highest  atomic  state  refers  to  the  highest  ionization  level  to  be  considered  by  the  pro¬ 
gram.  State  1  means  the  program  considers  only  the  neutral  atom;  state  2  means  the  program 
considers  the  neutral  atom  and  the  first  ionized  atom,  and  so  on.  Without  modification,  the 
program  considers  up  to  state  5,  four  ionizations. 

The  program  then  a.sks  for  the  clement  molecular  weight  (in  amu),  the  dissociation  and 
ionization  energies  of  the  honnmuclcar  diatomic  molecule  (in  cV).  If  the  diatomic  molecule 
docs  not  form,  enter  a  large,  negative  dissociation  energy,  say  uq  =  —50  eV,  and  a  large, 
postive  ionization  energy,  say  «/  =  50  eV.  This  convinces  the  program  the  “diatomic" 
molecule  very  readily  dissociates,  i.e.  does  not  form,  and  cannot  ionize  (since  it  does  not 
exist). 

ist:  Co] 

Atomic  mass  (amu)? 

^Paiidamonium 
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[  l.OOOE+00] 

MW:  C  l.OOOE+OO] 


Dissociation  energy  of  diatomic,  homonuclear  molecule  (eV)  ? 

[  -S.OOOE+01] 
uD:  C  -5.000E+01] 

Ionization  energy  of  diatomic,  homonuclear  molecule  (eV)  ? 

C  5.000E+01] 
ul:  [  5.000E+01] 

Next  the  program  asks  for  the  characteristic  temperatures  of  rotation  and  vibration  for 
the  neutral  and  singly  ionized  homonuclear  diatomic  molecule,  at  the  ground  electronic  state. 
These  can  be  found  in  Huber  &  Herzberg  fl5|  for  many  diatomic  molecules.  If  no  diatomic 
molecule  exists,  the  characteristic  temperatures  must  still  be  non-zero,  since  the  program 
divides  by  them  to  get  the  rotational  and  vibrational  partition  functions. 

Cheuracteristic  temperatures,  at  ground  electronicstate  (K)  ? 

Keutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
t  l.OOOE+OO  l.OCOE+00  l.OOOE+OO  l.OOOE+OO] 

The  program  then  summarizes  the  information  just  typed  in  about  the  element  and  asks 
if  it  meets  the  users  satisfaction.  If  not,  the  user  has  the  opportunity  to  re-enter  it;  if  so,  the 
user  can  append  the  new  clement  to  the  file  ABX.in,  so  the  element  will  be  “known"  by  the 
program  the  next  time  the  elctuent  is  u.secl. 

charT:  [  l.OOOE+OO  l.OOOE+OO  l.OOOE+OO  l.OOOE+OO] 

Summary : 

Highest  atomic  state  (1+highest  ionization  state) 

[0] 

Atomic  mass  (amu) 

[  l.OOOE+OO] 

Dissociation  energy  of  diatomic,  homonuclear  molecule  (eV) 

[  -5.000E+01] 

Ionization  energy  of  diatomic,  homonuclear  molecule  (eV) 

[  5.000E+01] 

Characteristic  temperatures,  at  ground  electronicstate  (K) 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
C  l.OOOE+OO  l.OOOE+OO  l.OOOE+OO  l.OOOE+OO] 

Is  everything  ok?  (l=ok,other=redo) 

Cl] 

Good.  Append  data  to  file  (ABX.in)?  (0*no,l=ye8) 

[0] 
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Next,  the  program  searches  through  the  file  ABXY.in  for  the  compound  formed  by  the 
two  chosen  elements.  File  ABXY.in  contair-s  energies  of  dissociation  and  ionization,  and 
characteristic  temperatures  for  rotation  and  vibration,  for  heteronuclear  diatomic  molecules. 
If  the  program  cannot  find  that  compound,  it  prompts  the  user  for  the  requisite  information. 

Will  not  append  data  to  file... 

Oops!  Couldn’t  find  CPn  in  file  (ABXY.in) 

Make  a  new  compound?  (I=ye8) 

Cl] 

Let’s  make  a  new  compound... 

Reaction  energy  of  diatomic,  heteronuclear  molecule  (eV)  ? 

C  -5.000E+01] 

RAB;  [  -5.000E+01] 

Ionization  energy  of  diatomic,  heteronuclear  molecule  (eV)  ? 

[  5.000E+01] 

lAB:  [  S.OOOE+01] 

Characteristic  temperatures,  at  ground  electronicstate  (K)  ? 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation.  Vibration 
[  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

As  for  the  homomiclcar  diatomic  molecule,  if  the  two  elements  form  no  compound,  the 
user  enters  a  zero  reaction  energy.  The  characteristic  temperatures  for  rotation  and  vibration 
must  be  non-zero,  though. 

The  program  then  summarizes  the  information,  asks  if  it  is  satisfactory,  and  offers  to 
append  the  new  information  to  the  file  ABXY.in. 

TAB:  [  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Reaction  energy  of  diatomic,  heteronuclear  molecule  (eV) 

C  -5.000E+01] 

Ionization  energy  of  diatomic,  heteronuclear  molecule  (eV) 

C  5.000E+01] 

Characteristic  temperatures,  at  ground  electronicstate  (K) 

Neutral:  Rotation,  Vibration;  Singly  Ionized:  Rotation,  Vibration 
C  l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00] 

Is  everything  ok?  (l=ok,other=redo) 
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Cl] 

Good.  Append  data  to  file  (ABXY.in)?  (0*no,l®ye8) 

[0] 


Finally,  the  program  has  information  it  needs  about  the  atomic  and  molecular  species, 
and  so  it  contiues  with  the  rest  of  the  user  input. 

Will  not  append  data  to  file... 

Using  C  ,Pn,CPn 

Stoichiometric  Proportions  of  C  ,Pn? 

C  O.OOOE+00.  l.OOOE+00] 

Atomic  Spectra  Input.  The  program  rcarls  the  atomic  spectra  data  in  from  fihis  on 
di.sk,  which  must  adhere  to  this  naming  convention: 

filename  =  (element). (charge) .in 

where  (element)  represents  for  the  element  symbol,  and  (charge)  represents  the  charge  of 
the  atom.  For  example,  if  the  user  wanted  to  examine  Pn  up  to  the  fifth  ionization  state,  he 
would  need  to  enter  data  into  four  files,  named  Pn_0.  in  through  Pn_4.  in. 

Each  file  contains  data  in  the  following  fashion;  the  first  line  contains  the  ionization 
energy  of  the  atom,  in  inverse  centimeters.  Subsequent  lines  contain  pairs  of  numbers:  the 
degeneracy  (dimensionless)  and  the  electronic  excitation  energy,  also  in  inverse  centimeters. 
The  electronic  excitation  energies  increase  up  to  the  ionization  energy  level.  The  last  line 
contains  two  negative  numbers,  to  flag  the  computer  to  stop  reading  degeneracy,  energy  level 
pairs.  The  user  can  find  degeneracies  and  electronic  cxitation  levels  for  atomic  species  in 
Moorc[14]. 

3.7  Benchmarking  the  Code 

Internal  Checks.  The  user  can  try  to  fool  the  code  by  entering  two  identical  elements 
with  non-zero  total  nuclear  fractious.  For  example,  say  the  user  enters  nitrogen  and  nitrogen. 

Using  N  .N  .NN 

Stoichiometric  Proportions  of  N  ,N 
[  5.000E-01,  5.000E-013 

The  program  does  not  recognize  that  these  two  elements  are  different.  Call  the  two 
elements  N  and  N’.  If  the  program  works  properly,  the  program  should  predict  the  same 
thermodynamic  properties  fer  elements  N  and  N’,  since  they  really  are  the  same.  Also,  the 
program  should  predict  equal  quantities  of  the  diatomic  molecules  Nj,  N’2,  and  NN’.  The 
output  file  NN.out,  partially  listed  here,  shows  that  the  program  predicts  correctly. 
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Partition  Function:  N 


atomic, 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

4.105E+00 

8.937E+00 

5.836E+00 

l.OOOE+00 

O.OOOE+00 

1.685E+01 

3.669E+01 

vib 

l.OOOE+00 

l.OOOE+00 

1 . OOOE+00 

l.OOOE+00 

O.OOOE+00 

2.087E+03 

2.159E+03 

rot 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

1.768E+00 

1.890E+00 

Partition  Function:  N 

atomic. 

atomic. 

atomic , 

atomic. 

atomic. 

diatomic. 

diatomic. 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

4.105E+00 

8.937E+00 

5.836E+00 

l.OOOE+00 

O.OOOE+00 

1.685E+01 

3.669E+01 

vib 

1 . OOOE+00 

1 . OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

2.087E+03 

2 . 159E+03 

rot 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

O.OOOE+00 

1.768E+00 

1.890E+00 

Partition  Function:  NN 

atomic. 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic, 

diatomic, 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

el 

1.685E+01 

3.669E+01 

vib 

2.408E+01 

1.037E+02 

rot 

2.779E+00 

3.175E+03 

Right-hand-sides  of  Saha  Equation 

atomic, 

atomic. 

atomic , 

atomic. 

atomic , 

diatomic. 

diatomic, 

neutral 

+1 

+3 

+4 

neutral 

+1 

N 

1.204E-09 

8 . 032E-23 

O.OOOE+00 

O.OOOE+00 

O.OOOE+00 

2.319E-03 

2.114E-10 

N 

1.204E-09 

8.032E-23 

0 . OOOE+00 

O.OOOE+00 

0. JOOE+00 

2.319E-03 

2.114E-10 

NN 

2.319E-03 

2.114E-10 

Nuclear  fractions. y 

atomic , 

atomic. 

atomic. 

atomic. 

atomic. 

diatomic. 

diatomic, 

neutral 

+1 

+2 

+3 

+4 

neutral 

+1 

N 

1.390E-02 

1.027E-06 

5.069E-24 

O.OOOE+00 

O.OOOE+00 

1.620E-01 

2.104E-06 

N 

1.390E-02 

1.027E-06 

5.069E-24 

O.OOOE+00 

O.OOOE+00 

1.620E-01 

2.104E-06 

NN 

1.620E-01 

2.104E-06 

yh.ye.yh+ye.beta  5.139E-01  8.367E-06  5,139E-01  1.628E-05 


Nitrogen  Density  Compare  the  density  vs.  temperature  profiles  of  nitrogen,  at  one 
atmosiihcre  pressure,  in  Figures  (2)  and  (3)[11,  p.l38]  They  follow  each  other  closely,  except 
Cambol  docs  not  consider  the  Nj  si)ccics. 

Air  Mole  Fraction  Compare  the  mole  fraction  vs.  temperature  profiles  of  air  (80% 
nitrogen,  20%  oxygen),  at  one  atmosphere  pressure,  in  Figures  (4)  and  (5)fU,  p.86]  Cambel 
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Figure  3;  Equilibrium  chemical  composition  of  nitrogen  plasma  at  atmospheric  pressure. 
(Cambel,  p.l38,  used  without  permission) 


31 


considers  species  that  this  program  does  not:  N2O,  N2O,  0“.  These  species  effect  the 
electron  density,  so  the  program  results  for  electron  density  do  not  match  well  with  Cambel, 
even  at  10,000  K.  The  dominant  species  at  10,000  K,  N,  O,  and  N2,  match  within  10%  or 
so.  However,  Cambel  keej>s  the  density  constant  at  the  sea-level  value,  while  the  program 
keeps  the  pressure  constant  at  one  atmosphere. 


4000  6000  0000  <  0000 

T«mo«r«tun  <K) 


Figure  4:  Mole  fraction  equilibrium  chemical  composition  of  air  at  one  atmosphere  pressure. 


Figure  5:  Mole  fraction  equilibrium  chemical  composition  of  air  at  sea-level  density.  (Cambel, 
p.86,  used  without  permission) 


Carbon  Dioxide  The  user  should  use  Selph’s  “Isp”  Code  for  chamber  temperatures  below 
6000  K.  Selph’s  code  considers  all  chemical  reactions  and  all  relavant  species.  However,  the 
JANNAF  tables  stop  at  6000  K,  so  at  temperatures  above  6000  K,  Selph  holds  the  heat 
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capacity  of  the  gas  constant  and  extrapolates.  The  author’s  program  calculates  theoretically 
the  equilibrium  state  function  using  the  Saha  Equation,  but  neglects  polyatomic  molecules. 
At  temperatures  much  above  6000  K,  this  assumption  becomes  minor,  as  almost  all  of  the 
molecules  have  dissociated. 

Table  (1)  compares  the  dominant  species  densities  predicted  by  Selph’s  and  the  author’s 
codes,  for  one  part  carbon  and  two  parts  oxygen,  at  6000  K  chamber  temperature  and 
one  atmosphere  chamber  pressure  pressure.  Note  that  the  dominant  species  are  close,  but 
discrepancies  arise  to  give  swound  10%  difference  (for  example,  atomic  carbon,  diatomic 
oxygen,  and  carbon  dioxide).  Also  note  that  at  6000  K,  there  Selph  predicts  less  than  1% 
CO2. 


Table  1:  Particle  density  comparision  between  Selph’s  and  author’s  codes  for  CO2,  at  6000 
K  and  one  atmosphere.  _ 


Selph,  “Isp” 
(m-’) 

Author,  “AB” 
(■»-’) 

0 

6.1E+23 

5.7E+23 

CO 

5.9E+23 

6.2E+23 

c 

l.OE+22 

2.3E+20 

02 

8.0E+20 

2.4E+22 

C02 

2.2E+20 

O.OE+00 

Thermodynamic  Functions  Table  2  compares  thermodynamic  properties  for  nitrogen 
species  N,  N+,  and  N2,  using  the  JANNAF  tables  and  the  author’s  code.  Results  match 
closely  for  atomic  species,  but  differ  slightly  for  diatomic  species. 


Table  2:  Thermodynamic  properties  comparison  between  JANNAF  tables  and  author’s  pro¬ 


gram,  for  some  nitrogen  species. 


Cp,T  =  3000 

Cp,T  =  6000 

As,  6000  K, 
A/i,  6000  K, 


Cp,T  =  3000 

Cp,T  =  6000 

As,  6000  K, 
A/i,  6000  K, 


Cp,r  =  3000 

Cp,T  =  6000 

As,  6000  K, 
A/i,  6000  K, 


K  (J/mol/K) 

K  (J/mol/K) 
3000  I<  (J/mol) 
3000  K  (kJ/mol) 


K  (J/mol/K) 

K  (J/mol/K) 
3000  K  (J/rao!) 
3000  K  (kJ/mol) 


K  (J/mol/K) 

K  (J/mol/K) 
3000  K  (J/mol) 
3000  K  (kJ/mol) 


JANNAF 


20.98 

25.54 

15.63 

68.42 


20.97 

22.39 

14.95 

64.95 


113.20 


Author 


20.96 

25.52 

15.62 

68.40 


20.96 

22.37 

14.93 

65.00 


37.74 

43.42 

28.38 

124.20 
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#  A  FORTRAN  Source  Code 
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A.l  ab.for 

PROGRAM  AB 


901 

c 


c 

c 


90 


c 


IMPLICIT  REAL*8  (a-h.l-z) 

CHARACTER  chA*2.chB»2.naine*4,TorP*l 

INTEGER  ist A , istB , seein . seemid . maxout , densmol . hovmany 

REALMS 


Apressl , pres82 . temp , tempi , temp2 . delt . 

AMWA . DA2 , IA2 . TA2 (4) . ionA (5) . guA (5 . 200 . 2) . 
ftMWB.DB2.IB2,TB2(4) .ionB(5) ,guB(5.200.2) , 

4Z0A(7,3) .ZA(7,3) ,ubarA(7) ,uubA(7) ,rh8A(7) ,yOA.yA(7) . 
tZOB (7 , 3) . ZB (7 . 3) , ubarB (7) . uubB (7) . rhaB (7) . y OB . yB (7 ) . 
tenthA(7) .entrA(7) .CpAC7) , 

4enthB(7) ,entrB(7) .CpB(7) , 

4RAB,IAB,TAB(4), 

itZ0AB(2,3)  ,ZAB(2.3)  ,ubarAB(2)  ,uiibAB(2)  ,rh8AB(2)  ,yAB(2) . 
tenthAB(7) ,eiitrAB(7) ,CpAB(7) ,enthe,Cpe,entre,ye,yh,beta, 
ftenthal.heatcp.entrop.rnf A(7) ,mfB(7) ,mf AB(2) , 
ftdenA(7) ,deiiB(7) ,denAB(2) ,dirA(7) ,dirB(7) ,dirAB(2) , 
tmdot , Pabs , Pjet .thrust , Isp 

DATA  tmpO  /2.9815E2/ 

INCLUDE  format, for 

FORMAT  (5x, ’in.loop’ ,3x, ’wyhX-yOX’ ,3x, 'mid. loop’ ,6x, 'wyh-l' , 
it3x,  'out. loop’  ,3x,  'residual') 


WRITEC*,*) 'Program:  AB.FOR,  ver.  1,3;  Author:  ’// 
ft’R.Nachtrieb,  Jan  93;' 

WRITE(*,*)'OLAC-PL/RKFE  Edwards  AFB  CA  93523-5000’ 

WRITEC*,*) 'internet  email:  nachtriebOuiuc.edu’ 

WRITEC*,*) 

WRITE(*.*) 

-startup  function  gets  temperatures,  pressures,  spectroscopic 
-data  from  user  and  files 
CALL  startup ( 

tpres8l,pre8s2.TorP, tempi, temp2,howmany. 

ft chA , ist A . MUA , DA2 , I A2 , TA2 , guA . ionA , yOA , 

ftchB , istB . MVB , DB2 , IB2 . TB2 , guB , ionB , yOB , 

ftname ,  RAB ,  lAB ,  TAB  .seein ,  seemid  .meixout ,  mdot ,  densmol) 

-get  partition  functions  at  STP,  for  later  use  in  entropy  difference 
CALL  Z1 ( tmpO , ist A , guA , TA2 , ionA , ZOA , ubar A , uubA) 

CALL  Z1  (tmpO ,  istB ,  guB , TB2 ,  ionB . ZOB ,  ubairB , uubB) 

CALL  Z2(tmp0 , ZOA , ubar A . uubA , ZOB , ub2u:B , uubb , TAB , ZOAB , ubar AB . uubAB) 
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c - either  search  for  P(T)=Pwant  or  increment  temperture 

IF  (TorP.EQ.'PO  THEM 
c - define  temperature  bounds 


Tmins6e3 

Tmax=30e3 

Tlo“Tmin 

Thi=Tmax 

resP=le7 

epsP=le-2 

c - define  lower  absorbed  power  bound 

WRITEC*,*) ’Getting  lower  absorbed-power  bound...’ 

WRITEC*,*) 

temp^Tmin 

CALL  saha( 

k  pressl . temp , ist A . istB , seein . seemid.mazout , 
k  HVA.DA2.IA2.TA2,ionA.guA. 
k  HWB.DB2.IB2.TB2.ionB.guB. 
k  ZOA.ZA.ubarA.uubA.rhsA.yOA.yA. 
k  Z0B.ZB.ubarB.uubB.rhs3.y0B.yB. 
k  enthA.entrA.CpA, 
k  enthB.entrB.CpB. 

k  RAB . I AB . TAB . ZOAB . ZAB . ubarAB . uubAB . rhs AB . yAB . 
k  enthAB.entrAB.CpAB.enthe.Cpe.entre.ye.yh.beta) 

CALL  gas( 
k  pressl. temp. 

k  DA2 . IA2 . ionA . yA . enthA , CpA . entrA . 
ft  DB2 . IB2 . ionB . yB . enthB . CpB . entrB . 

ft  name , RAB . I AB . yAB . enth AB , CpAB . entr AB . enthe . Cpe . entre . 
ft  enthal . heatcp . entrop . mf A , rof B . mf AB , denA , denB . denAB . 
ft  dirA.dirB.dirAB) 

Pabssenthal+mdot’t'leS/CyOA^MWA+yOB+MWB)  !W 
Pmin=Pabs 

c - define  lower  absorbed  power  bound 

WRITEC*.*) ’Getting  upper  absorbed-power  bound...’ 

WRITEC*.*) 

temp=Tmax 

CALL  sahaC 

ft  pressl. temp. istA. istB. seein. seemid.maxout. 
ft  MWA.DA2.IA2,TA2,ionA.guA. 
ft  MWB.DB2.IB2.TB2,ionB.guB, 
ft  ZOA.ZA.ubarA.uubA.rhsA.yOA.yA. 
ft  ZOB.ZB.ubarB.uubB.rhsB.yOB.yB. 
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ft  entliA,entrA,CpA. 
ft  entliB,entrB,CpB, 

ft  RAB . lAB . TAB . ZOAB . ZAB . ubar AB . uttbAB . rbsAB . yAB , 
ft  eatbAB , entr AB . CpAB , entbe , Cpe . entre . ye , yb , beta) 

CALL  gasC 
ft  presel.temp, 

ft  DA2 , IA2 , ionA , yA , enthA , Cp A , entr A . 
ft  DB2 , IB2 , iouB , yB , entbB , CpB , entrB , 

ft  name . RAB . TAB . yAB , entbAB , CpAB . entr AB , entbe . Cpe , entre , 
ft  entbal , beatcp , entrop . mf A ,mf B , mf AB , denA . denB , denAB . 
ft  dirA.dirB.dirAB) 

Pab8=enthal*mdot*le3/(yOA*MWA+yOB*MWB)  !W 

Pmaz=Pabs 

WRITECe.*) 

WRITE(*.*(A,1P2(E11.3))')»  Pmin.Pmax  (W)  » .Pmin.Pmax 
WRITE(*,*)'Wbat  power  to  absorb?  (W)  * 

READ(*,*)Pwant 

c - check  if  requested  power  within  power  bounds 

IF  (Pwant.LT.Pmin)  THEN 
WRITE (*,*) ’Sorry,  requested  power  too  low’ 
WRITEC*,*) ’Try  decreasing  the  mass  flow  rate’ 
WRITEC*,*) 

GOTO  90 

ELSEIF  (Pwant.GT.Pmax)  THEN 

WRITEC*,*) ’Sorry,  requested  power  too  high’ 
WRITEC*,*) ’Try  increasing  the  mass  flow  rate’ 
WRITEC*,*) 

GOTO  90 
ENDIF 

c - perform  temperature  search 

WRITEC*,*)’  resP  epsP  temp', 

ft  ’  out. loop  residual’ 

WRITEC*,*) 

DO  WHILE  CresP.GT.epsP) 

IF  CPabs.LT.Pwant)  THEN 
Tlo=temp 

t  emp=0 . 5*  C  t  emp+Th i ) 

ELSEIF  CPwant.LT.Pabs)  THEN 
Thi=temp 

temp=0 . 5*  CTlo+temp) 

ENDIF 
CALL  sahaC 
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993 


ft  pressl , temp , istA , istB , seein , seemid ,maxout , 

ft  MWA.DA2.IA2,TA2,ionA.guA. 

ft  MVB.DB2.IB2.TB2.ionB.guB. 

ft  ZOA.2U.ttbarA.uabA.rh8A.yOA.yA. 

ft  ZOB.ZB.ubarB.ttubB.rhsB.yOB.yB. 

ft  enthA.entrA.CpA. 

ft  enthB.entrB.CpB, 

ft  RAB . I AB . TAB . ZOAB , ZAB . ubar AB . uubAB . rhs AB . yAB . 

ft  enthAB , entr AB . CpAB . enthe . Cpe . entre . ye , yh , beta) 

CALL  gasC 
ft  pressl.temp. 

ft  DA2 . I A2 , ionA , yA , ent hA . CpA . entr A . 

ft  DB2 . IB2 . ionB . yB . enthB . CpB . entrB . 

ft  name . RAB . I AB . yAB , enthAB . CpAB . entr AB . enthe . Cpe . entre . 

ft  enthal .heatcp . entrop .mf A .mf B .mf AB .denA . denB . denAB . 

ft  dirA.dirB.dirAB) 

Pabs=enthal*mdot*le3/ (yOA*MWA+yOB*MWB)  ! W 
re3P=DABS((Pabs-Puant)/Pwant) 

F0RMAT(»+',1P3(E11.3)) 

WRITEC* , 993) resP . epsP . temp 
WRITE(*,*) 

ENDDO 

WRITEC*.*) 

WRITE(* . ' ( A , 1P3 (El 1.3))’)’  Pwant . Pabs . temp ’ , Pwant , Pabs , temp 

hovmany=l 

itemp=l 

— -calculate  ideal  rocket  performance 
CALL  rocket ( 

ft  yOA , yOB , HWA . MWB , ye , yh , temp , pressl , press2 ,mdot , 
ft  thrust, Isp.Pjet) 

CALL  output ( 

ft  yOA . yOB . tempi . temp2 . howmany . pressl . pressZ.mdot . temp . itemp , 
ft  chA . ZA . rhs A , yA . enthA , CpA , entr A , 
ft  chB . ZB . rhsB , yB , enthB . CpB . entrB . 

ft  name . ZAB , rhs AB , yAB . enthAB . CpAB . entr AB .enthe , Cpe . entre . 
ft  enthal , heatcp , entrop , mf A .mf B . mf AB . denA , denB . denAB . 
ft  dirA , dirB , dirAB , Pabs , Pj  et . thrust , Isp. densmol) 

WRITEC*,*) 

ELSEIF  (TorP.EQ. ’T’)  THEN 
delt=0. 

IF  (howmany. NE. 1)  delt=(temp2-templ)/DFL0AT(howmany-l) 
temp=templ 
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DO  100  itempsl.bovmany 

WRITE(*,192)’  tempCK),  preasKPa) * .temp.pressl 
imiTE(*.901) 

WRITEC*,*) 

c - ““get  equilibrium  molecule/atom/ion  fractions  using  pairtition  function 

CALL  sahaC 

A  pressl.temp.istA.istB.seein.seemid.meucout , 

A  MWA,DA2,IA2,TA2.ionA,guA, 

A  MVB.0B2.1B2.TB2,ionB.guB. 

A  ZOAfZA.ubarA.uubA.rhsA.yOA.yA. 

A  ZOB.ZB,ubarB,uubB,rhsB,yOB,yB. 

A  enthA,entrA,CpA, 

A  enthB,entrB,CpB, 

A  RAB .  I AB ,  T  AB .  ZOAB .  ZAB .  ubar  AB ,  uubAB ,  rhs  AB ,  yAB . 

A  enthAB , entr AB , CpAB , enthe , Cpe , entre , ye , yh , beta) 

c - calculate  ideal  rocket  performance 

CALL  rocket ( 

A  yOA , yOB , MWA , MWB , ye , yh , temp , pressl , press2 , mdot , 

A  thrust, Isp.Pjet) 

c - get  thermodynamic  properties  of  the  gas 

CALL  gasC 
A  pressl, temp, 

A  DA2 . I A2 , ionA , yA , enthA , CpA , entr A , 

A  DB2 , IB2 , ionB , yB , enthB , CpB , entrB , 

A  name , RAB , I AB , yAB , enthAB , CpAB , entrAB , enthe , Cpe , entre , 

A  enthal , heatcp . entrop , mf A ,mf B , mf AB , denA , denB , denAB , 

A  dirA,dirB,dirAB) 

c - calculate  maximum  absorbed  power 

Pabs=enthal*mdot*le3/(yOA*MWA+yOB*MWB)  !W 
c - ( J/mol)  (kg/s)  (g/kg)/Cg/mol)=W 

c - write  all  accumulated  data  to  file  for  later  digestion 

CALL  output ( 

A  yOA , yOB , tempi , temp2 , howmany , pressl , pressZ ,mdot , temp , itemp , 

A  chA , ZA , rhs A , yA , enthA , CpA , entrA , 

A  chB , ZB , rhsB , yB , enthB , CpB , entrB , 

A  name , ZAB , rhs AB , yAB , enthAB , CpAB , entrAB , enthe , Cpe , entre , 

A  enthal , heatcp . entrop ,mf A , mf B ,mf AB , denA , denB , denAB , 

A  dir A , dirB , dir AB , Pabs , P j  et , thrust , Isp , densmol ) 

WRITE(*.f) 
temp=temp+delt 
100  CONTINUE 

c - close  all  open  files,  tell  user  which  files  to  examine  for  output 
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ENDIF 

CALL  shutdown (hovnany , chA , chB , name . densmol ) 

END  imain  AB 

C^L*^cilt*^:^t1t4inf***nc*^^l****tlt************************************************ 

INCLUDE  STARTUP. FOR 
INCLUDE  GETX.FOR 
INCLUDE  GETXY.FOR 
INCLUDE  NEVX.FOR 
INCLUDE  NEWXY.FOR 
INCLUDE  FILENAME. FOR 
INCLUDE  OPUN.FOR 
INCLUDE  CLOZE. FOR 
INCLUDE  GETGU.FOR 
INCLUDE  21. FOR 
INCLUDE  22. FOR 
INCLUDE  SAHA. FOR 
INCLUDE  RHSl.FOR 
INCLUDE  RHS2.F0R 
INCLUDE  THERMl.FOR 
INCLUDE  THERM2.F0R 
INCLUDE  THERME.FOR 
INCLUDE  GESl.FOR 
INCLUDE  ITERY.FOR 
INCLUDE  GETYE.FOR 
INCLUDE  GETYH.FOR 
INCLUDE  BALANCE. FOR 
INCLUDE  NUY.FOR 
INCLUDE  GETWYH.FOR 
INCLUDE  RES 1. FOR 
INCLUDE  RES2.F0R 
INCLUDE  ROCKET. FOR 
INCLUDE  OUTPUT. FOR 
INCLUDE  GAS. FOR 
INCLUDE  SPECSHOW.FOR 
INCLUDE  ZISHOW.FOR 
INCLUDE  Z2SH0W.F0R 
INCLUDE  OPUNHDR.FOR 
INCLUDE  SHUTDOWN. FOR 
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A.2  startup.for 


Qi^t^tHHc^^^Hf^tnt************************************************************* 

SUBROUTINE  startupC 

ftpressl , pres82 , TorP , tempi , temp2 . howmany . 
ftchA . ist A . KVA , DA2 . I A2 . TA2 . guA , ionA , yOA . 
ftchB , istB , MWB , DB2 , IB2 , TB2 , guB , ionB , yOB , 

Aname , RAB . lAB . TAB . seein , seemid .maxout ,mdot . densmol) 

CHARACTER  chA*2 , chB*2 , chdum*2 , name*4 , TorP* 1 
INTEGER 
AistA, 

AistB, 

Aseein , seemid , maxout , densmol , howmany 
EEAL*8  press 1 , press2 , tempi , temp2 , 

AMWA ,  DA2 . I A2 . TA2 (4) . guA (5 , 200 . 2) , ionA (5) , yOA , 

AMWB, DB2,IB2.TB2(4).guB(S.200.2),ionB(5).y0B, 

ARAB,IAB,TAB(4) 

CHARACTER  str*30 
IMPLICIT  REAL*8  (a-h,l-z) 

INCLUDE  format. for 
UATA  atm  /0.101325e6/ 

c - Read  in  default  inputs  (stored  from  last  run;  if  no  file  AB.in  exists 

c - i.e.,  first  time  the  program  is  run,  program  goes  to  label  100  and 

c - starts  with  the  “hard-wired”  defaults. 

c - User  can  select  default  by  pressing  <Enter>,  or  enter  his  own  selection. 

c - After  entering  data,  choices  are  written  to  file  AB.in  for  next  time. 

OPEN  (UNIT=1.  FILE=' AB.IN') 

READ( 1 , * , ERR=100 , END=100) chA , chB 
READ ( 1 , * , ERR=100 , END=100) yOA , yOB 
READ( 1 , * , ERR=100 , END=100)pressl , press2 
READCl , * , ERR=100 , END=100)TorP 
RE AD ( 1 , * , ERR= 10C,END=100)templ,terap2, howmany 
READCl , ♦ , ERR=100 , END=100)den8mol 
READCl , * , ERR=100 , END=100) seein , seemid .maxout 
READCl . * . ERR=100 , END=100)mdot 
GOTO  200 
100  CONTINUE 
chA=’C' 
chB=’0' 
y0A=l . 
yOB=2. 
pressl=l. 
press2=0. 


44 


TorP=’T' 

templ-6e3 

temp2-30e3 

hovmanysl 

deo8mol=i 

seein=20 

seemids20 

]nazout=20 

mdot=l . 

200  CONTINUE 


c - get  element  choices;  case  sensitive  (thinks  lover  case  is  a  new 

c - element,  and  prompts  user  lor  element  data 


300  WRITE(*.*)*Ar,  C.  H.  N.  Q» 

WRITE(*,*) ’Which  two  (2)  elements  to  use?’ 
WRITE(*,’(”  C”.A,”,”A,”]  ”)’)chA,chB 
READ(*.’(A)’)str 

IF  (str.NE.’  ’)  READ<str,*)chA,chB 
IF  (chA.GT.chB)  THEN 
chdum=chA 
chA=chB 
chB~chdum 
ENDIF 

OPEN  <UNIT=10.  FILE* 'ABX. in ’.STATUS* ’old’) 


c - look  in  file  ABX.in  for  information  about  generic  unknown  element 

c - “X”,  where  X  is  letter  user  typed.  If  cannot  find  "X",  in  file 

c - ABX.in,  get  bomb!=0,  and  user  has  to  either  provide  information 

c - about  element  X,  or  select  new  elements  (start  over). 

c - When  data  is  satisfactorily  known  about  X,  data  is  assigned  to  (now) 

c - known  element  “A”,  which  the  rest  of  the  prograun  uses.  Process 

c - repeats  for  unknown  X  — >  known  B 


CALL  getX (chA , ist A . MWA . DA2 , IA2 , TA2 , bomb) 

IF  (bomb.NE.O.)  THEN 
iok*! 

WRITE(*,*) 

WRITE(*,*) ’Oops !  Couldn”t  find  ’,chA,’  in  file  (ABX.in)’ 
WRITE(*,*) ’Make  a  new  element?  (l=yes) ’ 

WRITE(*,’(”  C”.Il,”]”)’)iok 
READ(*,’(A)’)str 
IF  (str.NE.’  ’)  READ(str.*)iok 
IF  (iok.EQ.l)  THEN 


c - get  the  info  about  unknown  element  “X”.  When  finished,  user 

c - has  option  of  appending  new  element  data  to  file  ABX.in,  so 
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c 


- element  is  "known”  next  time  program  is  run. 

CALL  newXCchA , ist A ,MWA . DA2 , IA2 .TA2) 

ELSE 

REWIKD(UNIT=10) 

GOTO  300 
ENDIF 
ENDIF 

REWIND(UNIT=10) 

CALL  getX ( chB . is tB . HWB , DB2 . IB2 , TB2 , bomb) 

IF  (bomb.NE.O.)  THEN 
iok=l 

WRITEC*,*) 

WRITE(*,*) 'Oops!  Couldn”t  find  *,chB,’  in  file  (ABX.in)* 
WRITEC*.*) ’Make  a  new  element?  (l^yes)* 

WRITEC*, ’(”  [”.11,”]  ”)’)iok 
READC*, ’CA)’)str 
IF  Cstr.NE.'  ')  READC8tr.*)iok 
IF  Ciok.EQ.l)  THEN 

CALL  newX CchB , istB , MWB . DB2 , IB2 ,TB2) 

ELSE 

GOTO  300 
ENDIF 
ENDIF 

CL0SECUNIT=10) 

IF  CCchAC2:2).EQ. ’  ’) .AND. CchBC2:2) .EQ. *  ’))  THEN 
neune=chACl :  l)//chBCl :  1) 

ELSEIF  CchAC2:2).EQ. ’  ')  THEN 
name=chA  C 1 : 1 ) //chB 
ELSEIF  CchBC2;2).Eq.’  ’)  THEN 
name=chA//chBCl : 1) 

ELSE 

name=chA//chB 

ENDIF 

OPEN  CUNIT=10.  FILE=’ABXY.in’,STATOS=’old’) 

CALL  getXYCname.RAB.IAB, TAB, bomb) 

IF  Cbomb.NE.O.)  THEN 
iok=l 

WRITEC*,*) 

WRITEC*,*) ’Oops !  Couldn’’t  find  ’,name,’  in  file  CABXY.in)’ 
WRITEC*,*) ’Make  a  new  compound?  Cl=yes)’ 

WRITEC*, ’C”  [”,I1,”]  ”)’)iok 
READC*, ’CA)’)str 


46 


IF  (str.NE.*  ’)  READ(8tr,*)iok 
IF  (iok.EQ.l)  THEN 

CALL  newXY (name. RAB.IAB, TAB) 

ELSE 

GOTO  300 
ENDIF 
ENDIF 

CL0SE(UNIT=10) 

WRITEC*, ’(A) ') ’+Using  '//chA//» , *//chB//* . '//name 
WRITEC*,*) 

c - get  ratios  of  element  A:B  (parts  of  A  vs.  parts  of  B).  Converts 

c - to  fractions  of  nuclei. 

HRITE(*,*) 'Stoichiometric  Proportions  of  '//chA//' , '//chB//'?' 
WRITE(*.'("  C".1PE11.3,",".1PE11.3."]  ")')y0A,y0B 
READ(*.'(A)')str 

IF  (str.NE.'  ')  READ(str,*)yOA,yOB 

WRITE(*.192)'+Proportions  of  '//chA//' , '//chB.yOA.yOB 
IF((chA.EQ.chB).AND.(yOB.EQ.O.))  THEN 
y0A=0. 
yOB=l. 

ENDIF 

WRITE(*,*) 

c - Chamber  pressure  used  for  Saha  eqs.  Exit  pressure  used  for 

c - ideal  rocket  calculations. 

WRITE(*,*) 'Chamber,  exit  pressures  (atm)?' 

WRITE(*,'("  [",1PE11.3.".".1PE11.3,"]  ") ')pre8sl.press2 
READ(*,'(A)')str 

IF  (str.NE.'  ’)  READ(str,*)pre8sl,press2 
WRITE(*,192) '+pressl ,pres82  (atm) ' ,pressl,press2 
WRITE(*,*) 

c - User  specifies  temperature  (calculates  absorbed  power)  or 

c - users  specifies  absorbed  power  (searches  for  corresponding  temp) 

500  WRITE(*,*) 'Specify  Temperature  (T)  or  Absorbed  Power  (P)?  ' 
WRITE(*.'("  [",A,"]  ")')TorP 
READ(*,’(A)')str 
IF  (str.NE.'  ')  READ(str.*)TorP 
IF  ((TorP.EQ. 'p’).OR.(TorP.Eq.'P'))  THEN 
TorP='P> 

ELSEIF  ((TorP.EQ. 't’).0R. (TorP.EQ. 'T'))  THEN 
TorP='T’ 

ELSE 

WRITE(*,*)'What  the...  Try  again!' 
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GOTO  500 
ENDIF 

501  FORMATC’+'.A) 

WRITE(*.501)TorP 
IF  (TorP.EQ.'TO  THEN 

c - get  chamber  temperature  range  (independent  variable  iterated) 

c - If  howmany=l,  program  operates  on  first  tempature  only,  regardless 

c - of  upper  temperature. 

WRITE(*,*) ’Chamber  temperature  range: 
ft  ' tempi (K) ,temp2(K) .howmany?’ 

WRITE(*.'("  [”.1PE11.3.”,’MPE11.3.",*M5,"]  ”)’) 
ft  tempi, temp2,howmany 
READ(*,»(A)»)str 

IF  (str.NE.’  ’)  READ(8tr,*)templ,temp2,howmany 
HRITE(* . ' (A . 1P2 (IPEll . 3) . 15) ’ ) 
ft  ‘+templ(K) ,temp2(K) .howmany’ , tempi, temp2,howmany 
ENDIF 

c - all  interesting  information  is  printed  to  file  if  running  only 

c - one  temperature;  if  multiple  temperatures  are  run,  program  prints 

c - output  data  to  separate  files;  due  to  file  limitations  of  MS-DOS, 

c - user  has  to  choose  between  molefraction  and  density  outputs.  Sigh... 

IF  (howmany. GT.l)  THEN 
WRITE(*,f) 

WRITE(*,*) 'Print  density  (1)  or  mole  fraction  (2)  ?’ 

WRITE(*,’("  [".II,"]  ")’)densmol 
READ(*,’(A)')str 

IF  (str.NE.’  ’)  READ(str,*)densmol 
WRITE(*, ’("+densmol  [".II."]  ")')densmol 
EF'TF 
WRITE(*.*) 

c - On  personal  computers,  saha  equations  may  take  a  long  time  to  converge, 

c - especially  if  no  numeric  co-processor  is  used.  To  entertain  the  user, 

c - and  reassure  them  that  the  program  is  running  (instead  of  hanging). 

c - variables  seein.seeraid,  and  maxout  count  the  number  of  times  the  inner, 

c - middle,  and  outer  nested  while  loops  proceed  before  pooping  out  to  the 

c - screen  the  loop  count  and  the  convergence .  If  the  computer  is  fast, 

c - printing  out  to  the  screen  every  loop  really  slows  the  progr<un  down; 

c - but  if  the  computer  is  slow,  the  user  can't  tell  if  the  program  is 

c - working  or  not. 

c - Sugesstions;  fast  machines,  set  seein,seemid,maxout=100, 100,20 

c - slow  machines,  seein,seemid,maxout=20,20,20. 

c - The  program  is  "done"  when  the  final  residual  (sum  of  all  the 
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c - absolute  values  of  the  saha  equation  Ihs-rhs)  is  less  than  epsilon, 

c - where  epsilon=le-7.  When  the  outer  loop  residual  gets  down  around 

c - le-6,  it  starts  bouncing  around,  le-6  is  still  pretty  small,  so  if 

c - the  program  isn’t  done  before  maxout  outer  loops,  the  program  stops. 

WRITEC*,*) 'View  loop  progress:  inner, middle, outer?’ 

WRITE(*.*(”  [".14, ”,”,14, ”,”,14, ”3  ”)») 
ftseein.seemid, maxout 
READ(*.'(A)')str 

IF  Cstr.NE.’  ’)  READ(8tr,*)seein,8eemid, maxout 
WRITE(*,113) ’+seein, seemid, maxout' ,  8eein,8eemid, maxout 
WRITEC*.*) 

c - for  rocket  stuff 

WRITEC*, ♦) ’mass  flow  rate  (kg/s)?’ 

WRITEC*. ’(”  [’’.1PE11.3.”]  ”)’)mdot 

READ(*.’(A)’)str 

IF  Cstr.NE.’  ’)  READC8tr.*)mdot 

WRITEC*, 111) ’+mdot: ’ ,mdot 

WRITEC*,*) 

REWIND CUNIT=1) 

c - write  info  to  file  AB.in  for  next  time  defaults 

WRITE(l,*)chA,'  ’.chB 
WRITEC l,*)yOA,yOB 
WRITECl,*)pressl,press2 
WRITE(l,*)TorP 

WRITEC 1 ,*) tempi , temp2 , howmany 

WRITEC 1 , *)densmol 

WRITECl ,*)seein, seemid, maxout 

WRITECl,*)radot 

CL0SECUNIT=1) 

c - convert  atm’s  to  Pa’s 

press l=pressl*atm 

press2=press2*atm 

total=yOA+yOB 

yOA=yOA/total 

yOB=yOB/total 

c - open  appropriate  data  files  and  get  atomic  spectroscopic  data; 

c - load  spectroscopic  data  into  array  neune  "gu’’:  g  for  degeneracy, 

c - u  for  electronic  excitation  energy  level;  ionization  energies, 

c - in  cm‘-l,  go  into  array  "ion” 

CALL  opunCchA, istA) 

CALL  getguCistA,chA,guA,ionA) 

CALL  cloze 
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CALL  opun(chB,istB) 

CALL  getguCistB.chB.guB.ionB) 

CALL  cloze 

END  leubroutine  startup 

Qj^^^^iJtiim^^tJti*********************************************************** 
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A. 3  getx.for 

c — get  data  from  file  ABX.in  (unit  10)  for  unknown  element  X: 
c — using  element  name  (ch),  gets  highest  element  state  to  consider,  i.e. 
c — highest  ionization  level  +  1  (ist),  gets  molecular  weight  (HV.  in 
c — g/mol),  gets  disso'' j.tion  (uD)  and  ionization  (ul)  energies  of  the 
c — diatomic, homonuclear  molecule,  in  eV,  and  gets  the  characteristic 
c — temperatures  of  rotation  and  vibration  for  the  neutral  and  singly 
c — ionized  diatomic  molecule  (in  K).  If  cannot  find  ch  in  the  file, 
C”set8  bomb  to  non-zero,  which  flags  the  program  outside  this  routine. 
C********************************************************************** 
SUBROUTINE  getX( 
ftch ,  ist ,  MV ,  uD ,  ul ,  chaorT ,  bomb  ) 

CHARACTER*2  ch.readch 
INTEGER  ist 

REALMS  MW.uD,uI,charT(4),bomb 
bomb=0 . 

READdO,*) 

READ (10 , * . ERR=100 , END=100)readch, ist , MW. uD . ul , charT 
DO  WHILE  (readch.NE.ch) 

READdO,  *,ERR=100.END=100)readch,  ist, MW, uD,uI.  charT 
ENDDO 
RETURN 
100  bomb=l . 

END 


A. 4  getxy.for 

c — using  the  diatomic  molecule  name,  composed  of  the  compressed, 
c — concatenated  element  names,  in  alphabetical  order,  subroutine  looks 
c— in  file  ABXY.in  for  information  about  the  molecule.  Gets:  reaction 
c — (RAB)  and  ionization  (lAB)  energy  of  diatomic,  heteronuclear 
c — molecule,  and  the  characteristic  temperatures  for  rotation  and 
c— vibration  for  the  neutral  and  singly  ionized  molecule  (both  at  ground 
c — electronic  state).  If  cannot  find  info,  bomb!=l.  If  molecule  is 
c — pseudo-heteronuclear,  e.g.  it  looks  for  diatomic  nitrogen,  even 
c — though  N2  is  homonuclear,  homonuclear  info  is  returned. 

Qlti^^**:^^^^l*****************************r|i**********’l^********************** 

SUBROUTINE  getXY( 

Aname , RAB , lAB ,TAB , bomb) 

CHARACTER*4  name,readname 
REALMS  RAB, I AB, TAB (4), bomb 
bomb=0. 

READdO,*) 

READ ( 10 , ♦ . ERR=100 . END=100)readname , RAB , lAB . TAB 
DO  WHILE  Creadname.NE.name) 

READ (10 , * , ERR=100 . END=100)readname , RAB , lAB , TAB 
ENDDO 
RETURN 
100  bomb=l. 

END 
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A.5  newx.for 

c — Prumpts  user  for  informatiou  about  new  elements.  Asks  for 
c— informat ion  that  would  otherwise  be  found  in  file  ABX.in.  After 
c — completion,  gives  user  opportunity  to  append  new  element  info  to  file 
C"ABX.in,  so  next  time  the  information  will  be  found.  Gets  highest 
c — ionization  state  (ist),  molecular  weight  (MW),  dissocation  (uD)  and 
c — ionization  (ul)  energies  of  diatomic  homonuclear  molecule,  in  eV,  and 
c — characteristic  temperatures  of  rotation  and  vibration,  at  the  ground 
c-electronic  state,  of  the  neutral  and  singly  ionized  diatomic  molecule. 

SUBROUTINE  newXC 
ftch, ist, MW, uD,uI, chart) 

CHARACTER  ch*2,str*30 
INTEGER  ist 

REAL*8  MW,uD,uI,charT(4).iok 
10  WRITE(*,*) ’Let’ 's  make  a  new  element...’ 
ist=0 
MW=1. 
uD=-S0. 
ul=50, 
charT(l)=l. 
charT(2)=l, 
charT(3)=l. 
chart (4) =1. 

WRITE(*,*) 

WRITEC*,*) 'Highest  atomic  state  (1+highest  ionization  state)?  ’ 

WRITE(*.’(”  [”.11,”] '’)’)ist 

READ(*, ’(A)’)str 

IF  (str.NE. ’  ’)  READ(8tr,*)ist 

WRITE(*.’(”+ist;  [”  ,11.  ”]  ”)’)ist 

WRITEC*.*) 

WRITEC*,*) 'Atomic  mass  (amu)?  ’ 

WRITEC*.  ’(”  [”,1PE11.3,”]”)’)MW 

READ(*,'(A)’)str 

IF  (str.NE.’  ’)  READ(str,*)MW 

WRITE (*,  ’(”+MW:  [”  .  IPEll .3.  ”]  ”) ’)MW 

WRITE(*,*) 

WRITE(*,*) ’Dissociation  energy  of  diatomic,  homonuclear  ’// 

& ’molecule  (eV)  ?’ 

WRITE(+,  ’(”  [”.lPE11.3,”]”)’)uD 

READ(*.’(A)')str 

IF  (str.NE.’  ’)  READ(str,*)uD 
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WRITE(*,’("+uD:  [' '.IPEll  .3,  •  *]  » ')’)uD 
WRITEC*,*) 

VRITEC*.*) ’looization  energy  of  diatonic,  homonuclear  ’// 
t’ttolecule  (eV)  ?* 

WRITE(*,’("  [*MPE11.3,"]»')’)uI 
READ(*,’(A)')str 
IF  (str.HE.’  *)  READ(8tr.*)uI 
WRITE(*,*(’'+uI;  [*MPE11.3.»']")’)uI 
WRITE(*,*) 

VnilTEC*,*) ’Characteristic  temperatures,  at  ground  electronic*// 
ft ’state  (K)  ?* 

WRITEC*,*) ’Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
ft’Rotation,  Vibration’ 

WRITEC*.  ’  (  ”  [’  ’  .1P4(E11.3) .  ’  ’]  ”)  ’)charT 
READ(*,’(A)’)str 

IF  (str.NE.'  ’)  READ(8tr,*)charT 
WRITE(*.’(”+charT:  [”  .1P4(E11.3),  ”]  ”) ’)charT 
WRITEC*,*) 

WRITEC*,*) ’Summary: ’ 

WRITEC*,*) ’Highest  atonic  state  (1+highest  ionization  state)  ’ 
WRITE(*.’(”  C”.Il,”]”)’)ist 
WRITE(*,*) ’Atomic  mass  (amu)  ’ 

WRITE(*.’(”  C”.1PE11.3,”]”)’)MH 

WRITE(*,*) ’Dissociation  energy  of  diatomic,  homonuclear  ’// 
ft’molecule  (eV)  ’ 

WRITE(*,’(”  [".IPEU.S,”]  ”)’)uD 

WRITEC*,*) ’Ionization  energy  of  diatomic,  homonuclear  ’// 
ft ’molecule  (eV)  ' 

WRITEC*, ’(”  [”,lPEH.3,”]”)’)uI 

WRITEC*,*) ’Characteristic  temperatures,  at  ground  electronic’// 
ft ’state  CK)  ' 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  '// 
ft’Rotation,  Vibration’ 

WRITEC*, ’C"  [”.lP4CE11.3),”]”)’)charT 
iok*l 

WRITEC*,*) ’Is  everything  ok?  Cl=ok,other=redo) ’ 

WRITEC*, ’C”  [”.Il,”]”)’)iok 

READC*, ’ CA) ’)str 

IF  Cstr.NE.’  ’)  READCstr,*)iok 

IF  Ciok.NE.l)  GOTO  10 

iok=0 

WRITEC*,*) ’Good.  Append  data  to  file  CABX.in)?  C0=no, l^yes) ’ 
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WRITE(*,»(»*  C’Ml.”]  »’)’)iok 
READ(*,»(A)')8tr 
IF  (str.NE.’  »)  llEAD(8tr,*)iok 
IF  (iok.EQ.l)  THEN 

HRITE(>I‘,4‘) 'Appending  data  to  file  (ABX.in)...’ 
WRITECIO , ' ( A , 12 , 1P7 (Ell . 3) ) ' )  ch . ist , HW , uD , ul , charT 
ELSE 

WRITE(*,*) 'Will  not  append  data  to  file...' 

ENDIF 

END  ! subroutine  neuX 
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A.6  newxy.for 

c — Asks  user  for  inforDation  about  unknown  compound.  Gets  reaction 
c — (RAB)  and  ionization  (lAB)  energy  of  heteronuclear  diatomic  molecule, 
c — in  eV.  Also  gets  characteristic  temperatures  of  rotation  and 
c — vibration,  at  the  ground  electronic  state,  of  the  diatomic  neutral 
c — and  singly  ionized  molecule.  User  can  append  data  to  end  of  input 
C"file  ABXY.in  to  save  typing  data  in  again  next  time. 

QMf********,lf*^f4i*^^:^^f:lf*:tf************************************************** 

SUBROUTINE  newXYC 
tname,RAB,lAB,TAB) 

CHARACTER  name*4,str'(‘30 
REAL*8  RAB, I AB. TAB (4) 

WRITEC*,*) 

10  HRITEC*,*) *Let”s  make  a  new  compound...' 

RAB=-50. 

IAB=50. 

TAB(1)=1. 

TAB(2)=1. 

TAB(3)=1. 

TAB(4)=1. 

WRITE (*,♦) 

WRITE(*,*) 'Reaction  energy  of  diatomic,  heteronuclear  '// 
ft 'molecule  (eV)  ?’ 

WRITE(*,'("  C",1PE11.3,”]  ")»)RAB 

READ(*,'(A)')str 

IF  (str.NE.'  ’)  READ(str,*)RAB 

WRITEC*. ' ( ' '+RAB:  ["  , IPEll .3, ’ '] ' ') ')RAB 

WRITEC*,*) 

WRITEC*,*) 'Ionization  energy  of  diatomic,  heteronuclear  ’// 

4 'molecule  (eV)  ?’ 

WRITEC*, ’("  C”.1PE11.3,"]  ")»)IAB 

READC*.’CA)')str 

IF  Cstr.NE.’  ’)  READCstr,*)IAB 

WRITEC*. 'C’+IAB:  ["  ,  IPEll  .3,  "]  ") ’)IAB 

WRITEC*.*) 

WRITEC*,*) 'Characteristic  temperatures,  at  ground  electronic'// 
ft 'state  CK)  ?' 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
ft’Rotation,  Vibration’ 

WRITEC*. ’C”  C’MP4CE11.3).”]”)’)TAB 

READC*, ’CA)’)8tr 

IF  Cstr.NE.’  ’)  READCstr,*)TAB 
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WRITE(* . » (  ”  +TAB :  [ » * , 1P4(E11 . 3) , * » ] ‘ ’ )TAB 
iok*l 

WRITEC*,*) 

WRZTEC*,*) ’Reaction  energy  ot  diatomic,  beteronuclear  *// 

A ’molecule  (eV)  ’ 

WRITE(*.’(”  C”.1PE11.3.”]”)’)RAB 

WRITE(*,*) 'Ionization  energy  of  diatomic,  hetaronucllIT  '// 
('molecule  (eV)  ’ 

WRITE(*,’(”  t”.iPE11.3,”]”)’)IA8 
WRITEC’f,*) 'Characteristic  temperaturtl,  at  ground 
ft 'state  (K)  ’ 

WRITEC*,*) 'Neutral:  Rotation,  Vibration;  Singly  Ionized:  ’// 
(’Rotation,  Vibration’ 

WRITEC*. ’(”  C”.1P4(E11.3).”]”)')TAB 
WRITEC*,*) 'Is  everything  ok?  Cl=ok,other»redo) ' 

WRITEC*, ’C”  C”,Il,”3’')’)iok 

READC*,’CA)’)str 

IF  Cstr.NE.’  ’)  READC8tr,*)iok 

IF  Ciok.NE.l)  GOTO  10 

iok=0 

WRITEC*,*) 'Good.  Append  data  to  file  CABXY.in)?  C0=no,l=yes) ’ 

WRITEC*. ’C”  C”.I1."]  ”)’)iok 

READC*, ’CA)')str 

IF  Cstr.NE.’  ')  READCstr,*)iok 

IF  Ciok.EQ.l)  THEN 

WRITEC*,*) 'Appending  data  to  file  CABXY.in)...’ 
WRITEClO.’CA,lP6CE11.3))’)naae,RAB,IAB.TAB 
ELSE 

WRITEC*,*) 'Will  not  append  data  to  file...’ 

ENDIF 

END 

C’l'*4ii«i*it<4:4i****************4'**4<4‘*:t‘**>l"t<4i«*4i««*4>**********«*4t****4'********>|i 
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A.7  opun.for 

c — open  the  necessary  files  to  get  the  atomic  spectra  data.  The  naming 
c — convention  is  as  follows:  the  input  file  names  start  with  ch  (removed 
c — of  spaces),  have  an  underscore  (**.”),  and  then  an  number,  0-4, 
c — which  corresponds  to  the  charge  on  the  atom.  Charge  -*-4  corresponds 
c — to  i8t»5,  e.g.  C»V,  carbon  ionized  four  times. 

C******************************^*****^t**m****************************** 

SUBROUTINE  opun( 
hch.ist) 

CHARACTER  ch4>2,filenum*l,f ilename*12 
INTEGER  i 
DO  10  i=0,ist-l 
WRITE(filenttm,’(Il)»)i 
0PEN(UNIT=10+i, 

&  FILE=filename(ch,  ’_’//filenum//*  .inO, 
ft  STATUS=’0LD') 

10  CONTINUE 

END  '.subroutine  opun 

C*********^f*******************************^nlf*******if************4t^fiHiH^** 
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A.  8  cloze.for 


c — close  all  the  opened  files 

Qiti***********:lf******************^fif***JHt***1^*****m******1f*************** 

SUBROUTINE  cloze 
CL0SE(UNIT=10) 

CL0SE(UNIT=11) 

CL0SE(UMIT=12) 

CL0SE(UNIT=13) 

CL0SE(UNIT=14) 

END  'subroutine  opun 
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A.9  getgu.for 

c — necessary  files  are  already  opened,  so  reads  from  files  unlt=10 
c— through  unit=10+ist-l  to  get  the  ionization  energy  of  the  atom,  in 
c— cm‘‘-l,  and  the  (degeneracy,  energy  level)  pairs  (energy  levels  also 
c — in  cm“-l).  Ionization  energies  stored  in  array  “ion",  and 
c — degeneracy,  level  pairs  stored  in  array  “gu“.  For  example, 
c — guA(5,200,2)  stores  the  200x2  (g,u)  pairs  for  the  6  species  of 
c — element  A. 

Q:ti^cilHf:^**Mfitl************************************************************* 

SUBROUTINE  getgu( 

Aist,ch,gu,ion) 

INTEGER  ist 

CHARACTER  ch*2,f ilenum*l,filename*12 
REAL*8  gu(5,200.2),ion(S) 

INTEGER  i.j.k 
IMPLICIT  REAL*8  (a-h,l-z) 

DO  100  i=l,5 
DO  100  j=l,200 
DO  100  k=l,2 

100  gu(i,j,k)=0. 

DO  200  i=l,i8t 
WRITE(f ilenum, » (II) ' ) i-1 
WRITE(*,*) 'Getting  data  from  file  ’// 
ft  filename(ch, ’.’//filenum//’ .in’) 
idisk=9+i 

READ(idisk,*)  ion(i) 
j=0 

READddisk,*)  degen, englev 
DO  WHILE  (degen. GT.O.) 
j=j+l 

gu(i,j,l)=degen 
gu(i,j,2)=englev 
READ(idisk,*)  degen, englev 
ENDDO 

200  CONTINUE 

END  (subroutine  getgu 

C*:|<%***i|i*4<***«*4i*4<it‘****iiiii<**«***>)<*>ti4!4i4E****>l>4i4i*4i4i***itt4<*:4‘)4<V*******4<«***>t^>|:# 
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A.IO  saha.for 


C"The  meat  of  the  program.  Subroutine  has  3  nested  while  loops,  one 
c — of  which  is  hidden  in  SUBROUTINE  balance.  Outer  loop  updates 
c — electron  and  heavy  fractional  densities  (ye.  yh) .  and  beta 
c— (beta=ye/(yh+ye)) .  Middle  loop  holds  ye, yh, beta  constant,  but 
c — updates  heteronuclear  reaction  products  AB,  AB+,  Inner  loops  update 
c — atomic  and  homonucleax  diatomic  species.  (See  SUBROUTINE  balamce  for 
c — inner  loops.) 

C**************^^****************^^*^^***^^*>l^********1‘********************* 

SUBROUTINE  saha( 

ftpressi.temp,istA,istB.seein,seemid,maxout, 
tMWA . DA2 , IA2 , TA2 , ionA . guA . 
ftHWB . DB2 , IB2 . TB2 . ionB . guB , 

AZO A , ZA , ubar A , uub A , rhs A , yOA , yA , 

&Z0B , ZB , ubarB , uubB , rhsB , yOB , yB , 

AenthA , entrA , CpA , 
feenthB , entrB , CpB , 

ARAB , I AB , TAB , ZO AB , ZAB , ubar AB , uubAB , rhs AB , yAB , 
tenthAB , entrAB , CpAB , enthe , Cpe , entre , ye , yh , beta) 

IMPLICIT  REAL*8  (a-h.l-z) 

INTEGER  ist A , istB , seein , seemid , maxout 
REALMS  press 1, temp, 

AMHA.DA2,IA2,TA2(4) , ionA (5) .guA(5,200,2) , 

4MWB,DB2, IB2,TB2(4) , ionB(5) ,guB(5,200,2) , 

&Z0A(7,3) ,ZA(7,3) ,ubarA(7) ,uubA(7) ,rhsA(7) ,yOA,yA(7) , 

&Z0B (7 , 3) . ZB (7 , 3) , ubarB (7) . uubB(7) , rhsB (7) , yOB , yB (7) , 

&enthA(7) ,entrA(7) ,CpA(7) , 

&enthB(7) ,entrB(7) ,CpB(7) , 

ARAB, I AB, TAB (4), 

&Z0AB(2,3) .ZAB(2.3) .ubarAB(2) .uubAB(2) .rhsAB(2) ,yAB(2) , 

4enthAB(7) ,entrAB(7) .CpAB(7) , 

*tye,yh,beta 

C  local  variables  for  subroutine  saha 

INTEGER  i.imid,iout,ityA,ityB 
REAL* 8 

iepsilon, chi.yrain, 

&yeA,yhA,resA, 

&yeB,yhB,resB, 

4resAB,wyh, residual 

DATA  epsilon, chi, ymin  /l.E-7,0.5,l.E-26/ 

c - get  the  partition  functions  for  species  of  A,B,  and  AB  for  the 

c - given  temperature  and  pressure. 
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CALL  Zl(temp, i8tA,guA.TA2, ionA,ZA»ubarA,uubA) 

CALL  Z1 (temp . IstB , guB , TB2 , ionB , ZB , ubarB , uubB } 

CALL  Z2(temp , ZA , ubar A , uubA , ZB .ubarB , uubb , TAB , ZAB , ubar AB , uubAB) 

c - get  the  right-hand-sides  (rhs)  of  the  eaha  equations  (R_i  in  the 

c - report),  which  are  specified  for  a  given  temperature  and  pressure. 

CALL  rhs 1 (temp . press i , ist A . ionA . ZA , 0A2 , X A2 , HVA . rhs A) 

CALL  rhsl (temp , press 1 , i stB , ionB , ZB , DB2 , IB2 , MWB , r hsB) 

CALL  rhs2(temp , press 1 . ZA , ZB , ZAB , RAB . I AB . HWA . MWB , rhs AB) 

c - get  thermodynamic  properties  for  the  individual  species;  depend 

c - only  on  temp,  press,  partition  functions,  and  ubar,  uub;  can  get 

c - whole  gas  thermodynamic  properties  later  by  taking  the  weighted 

c - average  thermodynamic  properites,  where  the  weights  are  the  mole 

c - fractions . 

CALL  therml (temp , pressl , DA2 . I A2 , ionA . ZA , ZO A , ubarA , uubA , 
feenthA , CpA , entrA) 

CALL  therml (temp , pressl , DB2 , IB2 , ionB , ZB , ZOB , ubarB , uubB . 
ftenthB , CpB , entrB) 

CALL  therm2(temp, pressl, DA2,DB2, RAB, I AB, ZAB, 20AB, 
hubar AB , uubAB , enthAB , CpAB , entr AB ) 

CALL  therme(temp, pressl ,enthe,Cpe,entre) 

DO  200  i-i.7 
yA(i)=0. 
yB(i)=0. 

200  CONTINUE 

DO  300  i=l,2 
yAB(i)=0. 

300  CONTINUE 

c - guess  (“ges")  the  initial  values  of  the  atomic  and  homonuclear 

c - atomic  species  fractions  by  assuming  that  (1)  elements  A.B  can 

c - mix,  but  do  not  react  chemically,  (2)  a  given  element  has  two 

c - dominant  species  present  at  any  one  time,  all  others  being 

c - negligibly  (sp?)  small.  For  a  single  element,  or  for  mixtures, 

c - this  method  will  give  the  right  answer  to  about  5’/,,  and  zippy 

c - quick  too!  But  chemical  reactions  between  A  and  B  require  less 

c - elegant  methods... 

CALL  gesl(rhsA,y0A,yA) 

CALL  gesl (rhsB,y0B,yB) 

c - itery  finds  the  dominant  species  in  an  element,  and  sets  the 

c - variably  ity  to  that  species  position  in  the  element  array.  In 

c - subsequent  calculations,  operations  are  performed  using  that 

c - species  as  the  ‘‘row  pivot*'.  This  reduces  numerical  error, 

c - which  means  the  set  of  equations  converges  more  quickly  to  a  good 
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c - answer. 

CALL  itery(yA,ityA) 

CALL  itery(yB,ityB) 

c - functions  that  get  the  sum  of  the  atomic  fractions  for  the 

c - “heavy”  fraction  (yh),  and  the  weighted  average  atomic 

c - fractions  for  th^  electron  fraction  (ye),  where  the  weights  are 

c - the  species  charges.  At  this  point  in  the  progreun,  this  is  just 

c - our  first  guess  at  ye.yh.and  beta. 


CALL  getye(yA,yeA) 

CALL  getye(yB,yeB) 

CALL  getyhCyA.yhA) 

CALL  getyh(yB,yhB) 
ye=yeA+yeB+yAB (2) 
yh=yhA+yhB+yAB ( 1 ) +yAB (2) 
beta=ye/(yh+ye) 


c - calculate  the  maximum  possible  fraction  of  species  AB.  as  if  all 

c - of  one  species  or  another  combined  to  form  the  compound;  maximum 

c - is  modified  by  the  rhsAB,  which  sort  of  indicates  how  much  the 

c - species  AB  dissociates  as  the  temperature  increases. 

ymaxAB=0 . 5-DSQRT(0 . 25-yOA+yOB/(l .+rh8AB(l) ) ) 

iout=0 

residual^l. 

c - outer  while  loop;  finished  when  either  (1)  residuaKepsilon,  which 

c - means  the  y’s  we  have  satisfies  the  system  of  equations  nicely, 

c - or  (2)  the  loop  has  tried  maxout  times  to  converge,  and  hasn't 

c - done  it  yet,  and  so  probably  won’t,  maxout  is  definable  by  the 

c - user  in  the  input  section;  maxout=20  works  well,  and  reasonably 

c - balances  speed  vs.  efficiency.  This  loop  balances  charge. 


DO  WHILE  ((residual. GT. epsilon) .AND. (iout.LT. maxout)) 
yeO=ye 
yhO=yh 
imid=0 

c - reinitialize  the  upper  and  lower  bounds  for  yAB 

yhiAB=ymaxAB 
yloAB=0. 
wyh=-l . 


c - middle  while  loop:  is  done  when  the  weighted  sum  of  the  heavy 

c - atomic  fractions  (wyh)  is  within  epsilon  of  1,  where  the 

c - weights  are  the  number  of  nuclei  in  the  molecule  (2  for 

c - diatomic,  1  for  atomic).  This  is  essentially  the  conservation 

c - of  mass  equation.  Within  this  loop,  charge  is  held  temporarily 

c - constant. 
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DO  WHILE  (DABS(wyh-l.),GT. epsilon) 


yAB0=yAB(l) 

c - inner  loop(s).  one  each  for  the  elements  A  auid  B:  is  done  when 

c - the  weighted  sum  of  the  atomic  fractions  (sum  of  element 

c - nuclei)  is  within  epsilon  of  the  total  atomic  fraction,  yO; 

c - This  is  essentially  conservation  of  mass  for  each  element. 

CALL  balance (ye , yh , beta , epsilon , chi , ymin , yAB , seein , 
k  yOA.yA.ityA.rhsA) 

CALL  balance (ye , yh , beta, epsilon , chi , ymin , yAB , see in , 
k  yOB,yB.ityB,rhsB) 

c - done  with  inner  loops,  so  now  update  yAB,  yAB+, 

c - under-relaxing  using  bisection  method.  Converges  quickly  for 

(. - monotonic  functions. 


yAB(l)=yA(l)*yB(i)/(yh+ye)/rhsAB(l) 

yABl=yAB(l) 

IF  (yABl.GT.yABO)  THEN 
yloAB=yAB0 

yAB( 1 )= ( 1 . -chi) ♦yABO+chi*yhiAB 
ELSEIF  (yABl.LT.yABO)  THEM 
yhiAB=yABO 

yAB(l)=(l.-chi)*yABO+chi*yloAB 

ENDIF 

IF  (ye.EQ.O.)  THEN 
yAB(2)=0 

ELSEIF  (ye.GT.O.)  THEN 

yAB(2)=yAB(l)*rhsAB(2)/beta 

ENDIF 

c - get  new  weighted  sums  of  heavies  for  middle  loop  check 

CALL  getwyh(yA,yAB,wyhA) 

CALL  getwyh(yB,yAB,wyhB) 

wyh=wyhA+wyhB 

imid=imid+l 

922  FORMAT  (’+’ ,22x,lP2(E11.3)) 

c - the  progress-o-meter 

IF  (mod(imid,seemid) .EQ.O)  THEN 
WRITE(*,922)imid,wyh-l. 

ENDIF 

ENDDO 

c - get  new  electron,  heavy  fractions  for  the  outer  loop  check 

CALL  getye(yA,yeA) 

CALL  getye(yB,yeB) 

CALL  getyh(yA,yhA) 
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CALL  getyh(yB,yhB) 
ye=yeA-»-yeB+yAB  (2) 
yh=ybA+yhB+yAB ( 1 ) +yAB (2 ) 
yel®ye 
yhl=yh 

c - update  electron,  heavy  density  fractions  by  under-relajcing, 

c - using  bisection  method.  Again,  best  for  monotonic  functions 

c - (such  as  ue  have). 

ye=(l . -chi)»yeO+chi*yel 
yh=(l . -chi)*yhO+chi*yhl 
beta=ye/(yh+ye) 

c - get  partial  and  total  residuals  for  outer  loop  check 

CALL  res 1 (yOA , yA . rhsA , yAB , ye , yh , beta , res A) 

CALL  resl(yOB,yB,rh8B,yAB,ye,yh,beta,resB) 

CALL  res2 (yA , yB , yAB , rhs AB , ye , yh , beta , res A6) 

residual=resA+resB+resAB 

iout=iout+l 

944  FORMAT  ( '+' .44x, 1P2(E11 .3)) 

c - are  we  almost  there? 

WRITE(* , 944) iout , residual 
ENDDO 

END  (subroutine  saha 

Q^^^^.■|^*^^*i^^^1^1|i*1H^^^:|^1^1^1^t****^^**********************4t*********************^^* 
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A.ll  zl.for 


c — Make  the  partition  function  for  the  element,  given  the  atomic  spectra 
c — data  in  arrays  gu  and  ion.  and  the  temperature.  Also  calculates  the 
c — average  electronic  excitation  energy  level,  ubar,  and  the  average 
c—squared  energy  excitation  energy  level,  uub  (short  for  (u*u)bar) . 
c — These  are  useful  later  in  calculateing  thermodynamic  properties,  such 
c — as  heat  capacity.  Array  ZA(7,3),  holds,  for  exeunple,  partition 
c-function  data  for  the  7  states  (atomic  neutral,  atomic  +l,...+4, 
c — diatomic  neutral,  diatomic  +1)  of  element  A,  for  electronic, 
c — rotational,  and  vibrational  excitation  (Zrot  and  Zvib  =1  for  atomic 
c — species) .  Electronic  partition  functions  for  diatomic  species  are 
c—approximated  by  appropriately  averaged  atomic  partition  functions, 
c — For  example,  Zel,A2+  \approx  (ZelA)*(ZelA+) .  Remember,  partition 
c — functions  multiply,  not  add. 

C*********************************************!^*******ni******^***i^titi*im 

SUBROUTINE  Zl( 

fetemp , ist , gu , charT , ion , Z , ubar , uub) 

INTEGER  ist 

REAL+8  temp,gu(5,200.2).charT(4),ion(5),Z(7,3).ubar(7).uub(7) 
INTEGER  i.j 

IMPLICIT  REAL*8  (a-h,l-z) 
tempiav=l. 4387/temp 
DO  100  i*l,7 
DO  100  j=1.3 
100  Z(i,j)=0 

DO  200  i=l,ist 
summ=0 . 
sumu=0 . 
sumuu=0. 

j=l 

DO  WHILE  (gu(i,j,l).GT.0.) 
term  =  0. 

IF((gu(i, j ,2)*tempinv.LT.70. ) 
ft  .AND.(gu(i,j,2) .LT.ion(i)-temp/1.4387)) 
ft  term=gu ( i , j , 1 ) fDEXP ( -gu ( i , j , 2) *tempinv) 

summ  =suinm  +term 
surau  =sumu  +term*gu(i , j ,2) 

8umuu=sumuu+terin*gu(i,  j  ,2)*gu(i,  j  ,2) 
j=j+l 
ENDDO 

Z(i ,  l)=summ 
Z(i.2)=l. 


GG 


Z(i,3)-1. 
ubar ( i )=8uou/8UMm 
uub  (  i  )  s8Q]iiuu/8umffl 
200  COHTINOE 

Z(6.1)=Z(1,1)*Z(1.1) 

Z(6,2)=temp/charT(l) 

Z(6 , 3)=temp/charT(2) 
ubar(6)=ubar(l)+ubar(l) 

uub(6)=uttb(l)-ubar(l)*ubar(l)+uub(l)-ubar(l)*ubar(l) 

Z(7.1)=Z(l.l)*Z(2.1) 

Z (7 , 2)=temp/charT (3) 

Z(7,3)=teap/charT(4) 
ub2u:  (7  ) =ubar  ( 1) +ubar  (2) 

uub(7)=uub(l)-ubar(l)*ubar(l)+ttub(2)-ubar(2)*ubar(2) 
END  'subroutine  Z 


A.12  z2.for 


c — Gets  partition  functions  for  diatomic  species  (hence  the  2  in  Z2) . 
c — Approximates  the  electronic  partition  function  by  appropriately 
c — average  atomic  partition  functions.  For  example,  Zel,AB+  is  an 
c— average  between  Zel.A,  Zel.B-*-,  Zel.A-*-,  and  Zel.B.  Although  this  is 
c — not  rigorously  correct,  it  gets  in  the  ballpark;  I  don’t  worry  too 
c — much  about  diatomic  species  since  (1)  this  code  is  meant  to  be  run  at 
c — high  temperatures,  where  there  aren’t  many  diatomic  species,  and  (2) 
c — Selph’s  Isp  code  dof.s  a  great  (and  better,  too)  job  at  lower 
c — temperatures.  Also  gets  averaged  ubar  and  (u*u)bar  for  diatomic 
c — species.  Note  that  of  ZA=ZB,  ZAB=ZA2=ZB2,  as  should  be.  Thus, 
c — entering  in  two  identical  species  does  not  confuse  the  code. 

SUBROUTINE  Z2( 

Atemp , ZA , ubar A , uubA , ZB , ubarB , uubb ,TAB , ZAB , ubarAB , uubAB) 

REAL* 8  temp, 

*ZA(7.3) .ubarA(7) ,uubA(7) . 
ftZB(7,3) .ubarB(7) ,uubB(7) , 

&TAB(4) .ZAB(2,3) ,ubarAB(2) ,uubAB(2) 

IMPLICIT  REAL*8  (a-h.l-z) 

ZAB(1,1)=ZA(1,1)*ZB(1.1) 

ZAB(l,2)=temp/TAB(l) 

ZAB(1.3)=temp/TAB(2) 
ubarAB  '  1 ) =ubar A (l)+ubar8(l) 

uubAB ( 1 ) =uubA ( 1 ) -ubar A ( 1 ) *ubar A ( 1 )+uubB ( 1 ) -ubarB ( 1 ) *ubarB (1) 
ZAB(2.1)=0.5*(ZB(2,1)*ZA(1,1)+ZA(2,1)*ZB(1,1)) 

ZAB(2.2)=temp/TAB(3) 

ZAB(2,3)=temp/TAB(4) 

ubarAB (2) =0 . 5* (ubar A ( 1 ) +ubarA (2)+ubarB( 1 ) +ubarB (2) ) 
uubAB(2)=0.5*( 

ftuubA ( 1 ) -ubar A ( 1 ) *ubar A ( 1 ) +uubB ( 1) -ubarB ( 1 ) *ubarB ( 1 ) + 
ftuubA(2)-ubarA(2)*ubarA(2)+uubB(2)-ubarB(2)*ubarB(2)) 

END  .subroutine  Z2 
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A.  13  therml.for 


c — gets  tie  thermodynamic  functions  (enthalpy,  heat  capacity,  and 
c — entropy  difference  from  STP)  for  an  element’s  species,  and  sticks 
c — them  in  arrays  enth(7),  Cp(7),  and  entiv7).  Uses  the  element's 
c— partion  functions  at  STP  (ZO)  and  the  current  T,P  (Z).  See  Theory 
c — in  paper  for  derivation  of  thermodynamic  functions. 

Q:^li^J|f*^!**:|^}^^^f^^Jt^^t*4,^***^^*l|^,|l.^,^,^^******************************************** 

SUBROUTINE  thermK 

ktemp,pressl.uD,uI,ion,Z,ZO,ubar,uub,enth,Cp,entr) 

REAL*8  temp.pres8l,uD,ion(5),Z(7,3) ,Z0(7,3) ,ubax(7) ,uub(7) , 
tenth(7} ,Cp(7) ,entr(7) 

INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 

DATA  gasR, cminv , atm/8 . 31439 , 8 . 06549E3 , 0 . 101325e6/ 

DATA  tmpC  /2.981SE2/ 
tempinv=l . 4387/temp 
dision=0 . 5*uD*cminv 
DO  100  i=1.5 

IF  (Z(i.l).EQ.O)  THEN 
enth(i)=0. 

Cp(i)=0. 

entr(i)=0. 

ELSE 

enth(i)*gasR*temp*(2.5+(di8ion+ubar(i))*tempinv) 
Cp(i)=gasR*(2.5+(uub(i)-ubar(i)*ubar(i))*tempinv*tempinv) 
entr ( i ) =gasR* ( 2 . 5*DL0G (temp/ tmpC ) +ubar ( i ) *t  emp inv+ 
t  DLOG (Z ( i . 1 ) /ZO ( i . 1 ) ) -DLOG (pres  s 1/atm) ) 

dision=dision+ion(i) 

END  IF 
100  CONTINUE 

enth (6) =gasR*temp* (4 . 5+ubar (6) *tempinv) 
Cp(6)=gasR*(4.5+(uub(6)-ubar(6)*ubar(6))*tempinv»tempinv) 
entr (6) =gasR* (4 . 5*DL0G (temp/tmpO) +ubar (6) *tempinv+ 
tDLOG ( Z ( 6 , 1) /ZO (6 , 1 ) ) -DLOG ( pres  s 1 /atm) ) 
enth(7)=gasR*temp*(4. 5+(uI*cminv+ubar(7))*tempinv) 
Cp(7)=gasR*(4.5+(uub(7)-ubar(7)*ubar(7))*tempinv*tempinv) 
entr(7)=gasR*(4. 5*DL0G(te'np/tmp0)+ubar(7)*tempinv+ 
ftDLOG (Z(7,l)/Z0(7,l)) -DLOG (pre  s  s 1/atm) ) 

END 

Q*if:tfm*if****Hi**t******************************************************** 
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A. 14  therm2.for 


c—calculates  the  thermodynamic  functions  (h,c_p,8-s0)  for  the 
c — heteronuclear  diatomic  species.  Sticks  them  in  arrays  enthAB(2), 
c — CpAB(2),  and  entrAB(2). 

SUBROUTINE  thenn2( 

ftt  emp , press 1 , DA2 , DB2 , RAB , I AB , ZAB , ZOAB . 
ftubarAB , uubAB , enthAB , CpAB , entrAB) 

REAL*8  temp , press 1 , DA2 , DB2 , RAB , ZAB(2 , 3) , ZOAB (2,3), 
kubarAB(2) .uubAB (2) .enthAB (2) .CpAB (2) .entrAB (2) 

IMPLICIT  REAL*8  (a-h.l-z) 

DATA  gasR, cminv , atm/8 . 31439 . 8 . 06549E3 , 0 . 101325e6/ 

DATA  tmpO  /2.9815E2/ 
tempinv=l . 4387/temp 
dision=(0 . 5* (DA2+DB2) -RAB) *cminv 
enthAB(l)=gasR*temp*(4.5+(dision+ubarAB(l) )*tempinv) 
CpAB(l)=ga8R*(4.5+(uubAB(l)-ubarAB(l)*ubarAB(l))*tempinv*tempinv) 
entrAB  (l)=gasR*  (4. 5*DL0G(temp/tmp0)+abcLrAB(l)*tempinv+ 
tDLOGCZABCl . l)/Z0ABCl . l))-DLOG(pressl/atm) ) 
di8ion=(0.S=t‘(DA2+DB2)-RAB+IAB)*cminv 
enthAB(2)=gasR*temp*(4.5+(dision+ubarAB(2))*tempinv) 
CpAB(2)=gasR*(4. 5+(uubAB(2)”ubarAB(2)*ubarAB(2))*tempinv*tempinv) 
entrAB (2) =gasR* (4 . 5*DL0G ( t emp/tmpO)+ubarAB (2) ♦tempinv+ 

4DL0G(ZAB (2,1) /ZOAB (2,1)) -DLOG(pressl/atm) ) 

END 

C************************************^******************************t*ti 


A.15  therme.for 

c — get  thermodynamic  functions  (h,c_p,8“s0)  for  electrons. 

C********************************************************************** 

SUBROUTINE  thermeC 
Rtemp, pressl . enthe , Cpe . estre) 

REAL*8  t emp, pres si, enthe, Cpe, entre 
IMPLICIT  REAL+8  (a-h,l-z) 

DATA  gasR, atm/8 . 31439 . 0 . 101325e6/ 

DATA  tmpO  /2.9815E2/ 
enthe=2 . 5*gasR*temp 
Cpe=2 . 5*gasR 

entre=gasR* (2 . 5*DL00<temp/tmp0)-DL0G (press 1 /atm) ) 

END 
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A.16  rhsl.for 


c — gets  the  right-hand-sides  (rhs)  of  the  Saha  equation  for  use  in  the 
c — nuclear  fraction  (y)  calculations.  See  paper  theory  for  derivations 
c — Say  there  are  the  neutral  atomic  species,  4  ionized  atomic  species, 
c — and  the  tvo  diatomic  species:  7  species  per  element.  We  only  need  6 
c — equations,  and  therefore  6  rhs's,  to  solve  for  these  7  species, 
c — since  the  seventh  equation  is  conservation  of  mass.  Therefore,  for 
c — the  7  species  ve  have,  rhs(S)  will  be  left  empty  (0),  and  rhs (6)  and 
c — rh8(7)  correspond  to  the  diatomic  species  equations.  If  we  only  go 
c — up  to  A+3,  then  rhs (4)  and  rhs (5)  will  be  empty. 

C********************************************************************** 
SUBROUTINE  rhsK 

fttemp,pre8sl,ist,ion,Z,uD,uI,MW,rh8) 

INTEGER  ist 

REALMS  temp,pressl,ion(S) ,Z(7,3) ,uD,uI,HW,rhs(7) 

INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 

DATA  coefD.coefI,K_eV  /2. 59466E3, 3. 33381E-2. 11604.85/ 
tempinv=l .4387/temp 
tempre=temp*temp*DSQRT(temp)/pressl 
DO  100  i=l,7 
100  rhs(i)*0. 

DO  200  i=l,ist-l 
boltz=0. 

IF  (ion(i)*tempinv.LT.70.)  boltz=DEXP(-ion(i)*tempinv) 
rhs(i)=coef I*Z(i+l,l)/Z(i,l)*boltz*tempre 
200  CONTINUE 

mu=HW*MW/(MW+MW) 
boltz=0 . 

IF  (uD*K_eV/temp.LT.70.)  boltz=DEXP(-uD*K_eV/temp) 
rh8(6)=coefD*mu*DSQRT(mu)*boltz*tempre* 
ftZ(l,l)*Z(l,l)/ 

&Z(6,1)/Z(6,2)/Z(6,3) 

boltz=0. 

IF  (uI*K_eV/temp.LT.70.)  boltz=DEXP(-uI*K.eV/temp) 
rhs(7)=coef  I*boltz*tempre’t' 
ftZ(7,l)*Z(7,2)*Z(7,3)/ 
ftZ(6,l)/Z(6,2)/Z(6.3) 

END  ! subroutine  rhsl 

Qif^^tf*************  ***************************************  **************** 
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A.17  rhs2.for 


c — gets  right-hand-sides  of  Saha  equations  for  diatomic,  heteronuclear 
C”species.  See  paper  theory  for  derivation. 

SUBROUTINE  rhs2( 

fctemp , pressl , ZA , ZB . ZAB . RAB , I AB .HWA , MWB , rhs AB) 

REAL*8  temp , press 1 , ZA (7 . 3) , ZB (7 , 3) , ZAB (2 , 3) , RAB , I AB , 
ftMWA.MWB.rhsAB(2) 

IMPLICIT  REAL*8  (a-h.l-z) 

DATA  coefD.coefI.K_eV  /2.59466E3,3.33381E-2,11604.85/ 

tempre=temp*temp*DSQRT(temp)/pre8sl 

bolt2=0. 

IF  (RAB*K_eV/temp.LT.70.)  boltz=DEXPC-RAB*K_eV/temp) 
muAB=MWA*MWB/ (HWA+MWB) 

rhs AB (1) =coef D*muAB*DSQRT (muAB) ♦boltz*tempre* 
m(l.l)*ZB(l.l)/ 
mB<l.l)/ZAB(l,2)/ZABCi.3) 
boltz=0, 

IF  (IAB*K_eV/temp.LT.70.)  bolt2=DEXP(-IAB*K.eV/temp) 
rhsAB(2)=coef I*bolt2*tempre* 
feZAB (2 . 1) *ZAB(2 , 2) *ZAB (2 . 3) / 
ftZAB(l,l)/ZAB(1.2)/ZAB(l,3) 

END  ! subroutine  rhs2 
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A.18  gesl.for 

c — guesses  (hence  “ges*')  the  initial  species  distribution  for  single 
c — (hence  element,  based  on  the  assumption  that  only  two  species 

c — dominate  at  at  T.P.  This  assumption  linearizes  the  otherwise  very 
c — nonlinear:  equations,  and  gives  the  right  answer  to  within  5*/.  for  a 
c — single  element,  or  for  an  inert  (no  chemical  reactions)  mixture  of 
c — multiple  elements.  Lickity  split. 

c — The  conditions  for  checking  which  two  species  are  dominant  are 
c — these:  say  we’re  exeimining  N_+l,  N_+2,  and  N_+3  (abbrv.  Nl,  N2,  N3) 
c — The  switchover  point  from  considering  N1,N2  to  considering  N2,N3 
c — comes  when  there  are  equal  concentrations  of  Nl  and  N3,  Using  this, 
c — the  fact  than  Nl+N2+N3=N_heavy,  and  the  fact  that 
C"Nl+2N2+3N3=N_electrons  all  give  simple  relations  for  the  crossover 
c — points;  what  makes  the  crossover  points  unique  for  an  element  are 
c — the  values  of  the  right -hand-sides  of  the  Saha  equations,  which 
c — depend  on  the  molecular  weight,  partition  functions,  etc. 

Qi^**:(^>|^^^*M^*****^^**********************i‘********************************* 

SUBROUTINE  gesl( 
ftrhs.yO.y) 

REAL*8  rhs(7),y0,y(7) 

IMPLICIT  REAL+8  (a-h,l-z) 

IF  (rhs(6)*rhs(l) .LT.le-3)  THEN 
y(l)=yO*DSQRT(rhs(6)/(4.+rhs(6))) 
y(6)*y0-y(l) 

ELSEIF  (rhs(l)*rhs(2).LT.0.25)  THEN 
y(2)=yO/DSQRT(l+rhs(l)) 
y(l)=y0-y(2) 

ELSEIF  (rhs(2)*rhs(3) .LT.0.44)  THEN 

y(3)=y0*(-0.5+DSqRT(0.25+2.t‘rhs(2)/(l.+rhs(2)))) 

y(2)=y0-y(3) 

ELSEIF  (rhs(2)*rhs(3).GE.0.44)  THEN 
y(4)=yO*(-l.+DSQRT(l.+3.*rhs(3)/(l.+rhs(3)))) 
y(3)=y0-y(4) 

ENDIF 

CALL  getye(y,ye) 

IF  (ye.EQ.O.)  THEN 

y(2)=y(l)*rhs(l)+DSQRT(y(l)*y(l)*rhs(l)*rhs(l)+ 
k  y(l)*rhs(l)*(y(l)+y(6))) 

y(7)=y(6)*rhs(7)+DSQRT(y(6)+y(6)*rhs(7)*rhs(7)+ 
ft  y(6)*rh3(7)*(y(6)+y(l))) 

ENDIF 

END  ! subroutine  gesl 


Qit:^^^tiii^**i^t**if*******************************************<l:**************‘* 
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A.  19  itery.for 

c — find  which  nuclear  fraction  is  largest,  and  set  the  equations  to 
c — “pivot”  around  that  species. 

C**4i*****4i*****«:*<*>ti**:4<4:*«****>4i*********4>4i****’)‘*4‘******4>4<**«*4>******4<«*4> 

SUBROUTINE  iteryC 
iy.ity) 

INTEGER  ity 
REAL*8  y(7) 

IMPLICIT  REALMS  (a-h,l-z) 
ity=6 

IF  (y(l).GT.y(6))  THEN 
ity=l 

ELSEIF  (y(2).GT.y(l))  THEN 
ity=2 

ELSEIF  (y(3).GT.y(2))  THEN 
ity=3 

ELSEIF  (yC4).GT.y(3))  THEN 
ity=4 
END  IF 

END  ! subroutine  itery 

C**,H1r*******iHf*:^t*1Hci^:^t*********************t*************************** 
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A. 20  getyh.for 

c— get  the  SUB  of  the  nuclear  fraction8=the  "heavy”  fraction,  yh. 

SUBROUTINE  getyhC 
ty.yh) 

REAL*8  y(7).yh 
INTEGER  i 

IMPLICIT  REALMS  (a-h.l-z) 
yh=0. 

DO  100  i=1.7 
100  yh=yh+y(i) 

END  ! subroutine  getyh 


A.21  getye.for 

c — get  the  electron  fraction,  ye,  based  on  the  charge  weighted  average 
c — fraction 

SUBROUTINE  getyeC 
fty.ye) 

REAL+8  y(7),ye 
INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 
ye=0. 

DO  100  i-1.5 

100  ye=ye+(i-l)*y(i) 

ye=ye+y(7) 

END  ! subroutine  getye 

C**************************************************^i*<Hf*^****^*^-^i**if*** 
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A.  22  feetwyh.for 

c— get  the  mass-weighted  sum  of  the  heavy  nuclear  fractions.  Used  for 
c — conservation  of  mass  equations. 

Q********************************************************************** 

SUBROUTINE  getwyh( 
ftyiyAB.wyh) 

REALMS  y(7),yAB(2),wyh 
INTEGER  i 

IMPLICIT  REALMS  (a-h,l-z) 
wyh=0. 

DO  100  i=1.5 
100  wyh=wyh+y(i) 

wyh=wyh+2 .♦y(6)+2.*y(7) +yAB ( 1 ) ♦yAB (2 ) 

END  ! subroutine  getwyh 

Q****i^***Mf****^t*********7>i:i^^t******************************************** 
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A. 23  balance.for 


c — the  inner  loop.  Finishes  when  either  (1)  the  mass  weighted 
c — nucleau:  fraction  sum  (wyh)  is  within  epsilon  of  the  total  nuclear 
c — fraction,  yO,  or  (2)  the  sum  of  the  heavy  fractions  (yhh,  a 
c — different  name,  so  we  don't  screw  up  the  outer  loop)  is  very  very 
c — small.  This  protects  us  against  the  case  where  all  of  an  element 
c — combines  with  the  other  to  form  the  compound  AB;  otherwise,  the  routine 
c — will  keep  trying  to  reduce  the  p2u:tial  residual  by  making  the 
c — concentrations  smaller  and  smaller,  until  finally  we  get  an 
c — underflow  error. 

C*)^!******************************************************************** 

SUBROUTINE  balance ( 

&ye . yh , beta, eps ilon , chi , ymin , yAB , seein , yO , y , ity , rhs ) 

INCLUDE  format. for 
INTEGER  ity, seein 

REALMS  ye, yh, beta, epsilon, chi, ymin, yAB(2) ,y0,y(7) ,rhs(7) 

INTEGER  iin 

IMPLICIT  REAL+8  (a-h,l-z) 

loy=0. 

hiy=yO 

iin=0 

wyh=-l . 


yhh=l . 

c - the  inner  loop. 

DO  WHILE  ((DABS(wyh-yO) .GT. epsilon) .AND. (yhh. GT. ymin)) 

c - get  the  new  (“nu'')  fractions  (“y") 

CALL  nuy ( ity , yO , y , ye , yh , beta , rhs) 

c - then  get  the  new  weighted  heavy  fraction 

CALL  getwyh(y,yAB,wyh) 
c - under-relar/bisection  the  pivot  y 


IF  (wyh.GT.yO)  THEN 
hiysyCity) 

yCity)=(l .-chi)*y(ity)+chi*loy 
ELSEIF  (wyh.LT.yO)  THEN 
loy=y(ity) 

y(ity)=(l.-chi)*y(ity)+chi*hiy 

ENDIF 


c - get  the  new  weighted  heavy  fraction  with  the  updated  pivot  y 

CALL  getwyh(y,yAB,wyh) 

c - and  get  the  test-case  heavy  fraction 

CALL  getyh(y,yhh) 
iin=iin+l 
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c 


hov'8  the  progress? 

IF  (mod(iin,seein) .EQ.O)  THEN 
WRITE('«‘.192)  ’+» ,iin.wyh-yO 
ENDIF 
ENDDO 
END 

C*^i*****t************************************************************** 
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A.24  nuy.for 

c — given  the  pivot  y  (ity  is  the  array  index  for  the  y  array),  the 
c — total  nuclear  fraction  yO,  the  temporary  ye, yh, beta,  and  the  rhs 
c — array,  we  update  the  other  nuclear  fractions  for  an  element. 

Q  %  :tE  :|i  4l  :|i  :4i  %  4c  #  #  1)1  1):  *  *  *  *  41  *  *  «  *  *  *  *  *  % ’^  41  *  41 ’ll  4^  *  ■I' * 

SUBROUTINE  nuy( 
ftity,y0.y,ye,yh,beta,rhs) 

INTEGER  ity 

REAL*8  y0,y(7) ,ye,yh,beta,rhs(7) 

IMPLICIT  REALMS  (a-h.l-x) 

c - check  if  electron  density  is  zero,  and  should  still  be  zero.  Need 

c - to  check,  because  otherwise,  will  try  to  divide  by  beta,  ;hich  is 

c - zero  if  ye=0. 

IF  (ye.EQ.O.)  THEN 
y(2)=0. 
y(3)=0. 
y(4)=0. 
y(6)=0. 
y(7)=0. 

c - calculate  the  other  non-charged  species  from  the  pivot  species 

IF  (ity.EQ.6)  THEN 
y(l)=DSQRT(rhs(6)*y(6)*(yO-y(6))) 

ELSEIF  (ity.EQ.l)  THEN 

y(6)=2.*y(l)*y(l)/rhs(6)/(y0+y(l)) 

ENDIF 

c - find  the  other  nuclear  fraction  y’s  from  the  pivot  y.  The 

c - equations  are  only  slightly  non-linear  when  we  (temporarily)  hold 

c - the  electron  fraction,  ye,  and  beta  (as  well  as  yh)  constant. 

ELSEIF  (ye.GT.O.)  THEN 
IF  (ity.EQ.6)  THEN 

y(l)=DSqRT((yh+ye)*rhs(6)*y(6)) 

y(2)=rhs(l)/beta*y(l) 

y(3)=rhs(2)/beta*y(2) 

y(4)=rhs(3)/beta*y(3) 

y(5)=rhs(4)/beta*y (4) 

y(7)=rhs(7)/beta*y(6) 

ELSEIF  (ity.EQ.l)  THEN 
y(2)=y(l)*rhs(l)/beta 
y(3)=y(2)*rhs(2)/beta 
y(4)=y(3)*rhs(3)/beta 
y(5)=y(4)*rhs(4)/beta 
y(6)=y(l)4y(l)/rhs(6)/(ye+yh) 
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y(7)=y(6)*rhs(7)/beta 
ELSEIF  (ity.EQ.2)  THEN 
y(l)ssy(2)*beta/rh8(l) 
y(3)=y(2)*rb8C2)/beta 
y(4)=y(3)*rhs(3)/beta 
y(5)=y(4)*rhs<4)/beta 
y(6)=y(l)*y(l)/rh8(6)/(ye+yh) 
y(7)*y(6)*rh8(7)/beta 
ELSEIF  (ity.Eg.3)  THEM 
y(2)«y(3)*beta/rh8(2) 
y(l)=y(2)*beta/rb8(l) 
y(4)=y(3)+rh8(3)/beta 
y (5)=y (4) *rh8 (4) /beta 
y(6)=y(l)*y(l)/rh8(6)/(ye+yh) 
y(7)=y(6)*rhs(7)/beta 
ELSEIF  (ity.EQ.4)  THEN 
y(3)=y(4)*beta/rh8(3) 
y(2)=y(3)*beta/rhs(2) 
y(l)=y(2)*beta/rh8(l) 
y(5)=y(4)*rhs(4)/beta 
y(6)=y(l)*y(l)/rh8(6)/(ye+yh) 
y(7)=y(6)*rhs(7)/beta 
ENDIF 
ENDIF 

END  (subroutine  nuy 

^♦♦♦♦♦♦♦♦♦♦♦♦♦♦iti*********************************************  ********** 
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A.25  resl.for 


c — get  the  residual  for  an  element,  subroutine  can  be  used  for  either 
c — element;  add  residuals  from  both  elements,  and  hetero-diatomics  to 
c — get  total  residual. 

Q^^tc*^!*************************************************************^!**** 

SUBROUTINE  resl( 

tyO , y , rhs , yAB , ye , yh , beta , re  s idual) 

REAL-xB  yO , y (7) ,  rhs (7 ) ,  yAB C2) , ye , yh , beta ,  residual 
INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 

REAL*8  res (7) 

CALL  getHyh(y,yAB.wyh) 
res(i)=  wyh-yO 

res(2)=-rhs(l)*y(l)+beta*y(2) 

res(3)=-rhs(2)*y(2)+beta*y(3) 

res(4)=-rhs(3)*y(3)+beta*y(4) 

res(5)=-rhs(4)*y(4)+beta*y(5) 

res(6)=  y(l)*y(l)/(yh+ye)-rhs(6)*y(6) 

res(7)=-rhs(7)*y(6)+beta*y(7) 

residual=0. 

DO  100  i=1.7 

residual=residual+DABS(res(i)) 

100  CONTINUE 

END  ! subroutine  resl 

C**************>^***********************************i|^>lf^i^il^t^:*^t************ 
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A.26  res2.for 


c — get  residual  for  heteronuclear  diatomic  species. 

Q***^,*if*^t**^t*7tt*^l******Tll************************************************ 

SUBROUTINE  res2( 

tyk , yB , yAB , rhsAB , ye , yb . beta , residual) 

REALMS  yA(7) , yB(7) , yAB (2) , rh8AB(2} , ye , yh , beta, residual 
INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 

REAL*8  res (2) 

res(l)*  yA(l)*yB(l)/(yh+ye)-rh8AB(l)*yAB(l) 
res (2) =-rh8AB (2) *yAB ( 1 ) +beta*yAB(2) 
residual=0. 

DO  100  i=1.2 

residual=residual-t-DABS(res(i)) 

100  CONTINUE 

END  ! subroutine  res2 

Q***it,^t^t1ti):^t************************************************************ 


85 


A.  27  gas.for 

Qili^c**yl/i*1Hf)tf*****W:******************************************************* 

SUBROUTINE  gas( 

£pressl,temp, 

ftOA2 . IA2 , ionA . yA , enthA , CpA . entrA . 

ADB2 , IB2 , ionB , yB , enthB . CpB . entrB , 

tname , RAB , lAB , yAB , enthAB . CpAB . entrAB , enthe , Cpe , entre , 
ftenthal ,  heat  cp ,  entrop ,  mf  A ,  mf  B ,  sif  AB ,  denA ,  denB ,  den  AB , 

AdirA , dirB . dir AB) 

CHARACTER  n2une*4 
REAL*8 

&pre8sl.temp, 

ftDA2 , IA2 . ionA (5 ) . yA  C7 ) . enthA (7 ) . CpA (7) , entr A ( 7 ) . 
tDB2. IB2, ionB(5) .yB(7) .enthB(7) ,CpB(7) ,entrB(7) , 

ARAB. IAB,yAB(2) ,enthAB(2) .CpAB(2) ,entrAB(2) . 

Aenthe , Cpe , entre , 

Aenthal, heatcp, entrop, mfA(7)  ,jnfB(7)  ,mf  AB(2)  , 
AdenA(7),denB(7).denAB(2),dirA(7),dirB(7).dirAB(2) 

INTEGER  i 

IMPLICIT  REAL*8  (a-h.l-z) 

DATA  bolt2k/1.38066E-23/ 

DATA  JmoleV , Jmolcm/96484 .52,11. 96263/ 

c - JmoleV  (96484.52  J/mol/eV)  converts  eV  to  J/mol 

c - re-evaluate  the  heavy  and  electron  fractions,  just  t'  ’  safe. 

CALL  getye(yA,yeA) 

CALL  getye<yB,yeB) 

CALL  getyh(yA,yhA) 

CALL  getyh(yB,yhB) 

ye=yeA+yeB+yAB (2) 

yh=yhA+7hB+y AB ( 1 ) +yAB ( 2 ) 

pkt=press 1/ (boltzk*temp* (ye+yh) ) 

enthal=0. 

heatcp^O. 

entrop=0. 

c - calculate  the  heats  of  formation  for  the  species; 

c - dissociation  (“d"),  ionization  (“i”),  and  reactions  (“r") 

c - also  calculate  the  gas  thermodynamic  properties,  based  on  the 

c - species  properties  and  the  mole  fractions. 

dirA(l)=0.5*DA2* JmoleV  ! J/mol 
dirB(l)=0.5*DB2* JmoleV 
DO  100  i=1.5 

enthal=enthal+yA ( i ) *enthA ( i) +yB( i) *enthB ( i ) 
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heat  cp=heat  cp+yA ( i ) *CpA ( i )+yB ( i ) ♦CpB ( i ) 
entrop=entrop+yA ( i) ♦entr A ( i) +yB ( i ) *entrB ( i) 
dirACi-t-D^dirACD+ionACDoJmolcD  !  J/mol 
dirB(i+l)=dirB(i)+ionB(i)»Jmolcm 
mf A(i)=yA(i)/(yh+ye)  Imole  fraction 
mfB(i)«yB(i)/(yh+ye) 
denA ( i ) =yA ( 1) «pkt 
denB  ( i  )  *y  B  ( i  )  ♦pkt 
100  COMTIHUE 

DO  200  i*l,2 
enthal=enthal+ 

ft  yAB(i)*enthAB(i)+yA(i+S)*enthA(i+5)+yB(i+5)*enthB(i+S) 
heatcp=heatcp-<- 

ft  y AB ( i ) • CpAB (i)+yA(i+5)*CpA(i+5) +yB ( i+5 ) ♦CpB ( i+5 ) 
entrop=entrop+ 

ft  yAB  Ci)*entrAB(i) +yA ( i+5 ) *entr A ( i+5) +yB ( i+5) *ent rB ( i+5) 
mf A ( i+5) =y A ( i+5 ) / (yh+ye ) 
mf B ( i+5 ) =yB ( i+5 ) / (yh+ye ) 
denA ( i+5) =yA ( i+5) *pkt 
denB(i+5)=yB(i+5)*pkt 
mfAB(i)=yAB(i)/(yh+ye) 
denAB ( i ) =yAB ( i ) ♦pkt 
200  CONTINUE 
dirA(6)=0. 
dirB(6)=0. 
dirA(7)=IA2*JmoleV 
dirB(7)=IB2+JmoleV 
dir AB ( 1 ) = ( 0 . 5* ( DA2+DB2 ) -RAB) ♦ JmoleV 
dirAB(2)=(0.5*(DA2+DB2)-RAB+IAB)*JmoleV 
eiithal=enthal+ye*enthe 
heatcp=heatcp+ye*Cpe 
entrop=entrop+ye*entre 
enthal=enthal/ (ye+yh) 
heatcp=heatcp/(ye+yh) 
entrop=entrop/ (ye+yh) 

c  WRITE(*,'(A,1PE11.3)’)’  Enthalpies  (J/mol) '// 

c  tname//'  gas’.enthal 

c  WRITE(*, ’(A,1PE11.3)’)’  Heat  capacities  (c_p)  (J/K/mol)’// 

c  tname//’  gas’.heatcp 

END  ! subroutine  output 

C^i********************************************************************* 
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A.28  rocket. for 


c — get  performance  from  an  ideal  rocket,  assuming  (1)  our 
c — temperature  and  pressure  are  that  of  the  rocket  cheunber,  (2)  frozen 
c — flow  (no  change  in  equilibrium  species  composition  downstream),  (3) 
c — isentropic  (reversible)  expansion  of  gases,  (4)  haird  sphere  ideal  gas 
c — (so  forget  all  the  complicated  thermodynamic  functions  we  calculated; 
c — just  use  c_p=2.5R).  Frozen  flow  assumption  relieves  us  from  having 
c — to  calculate  coronal  equilbrium  at  exit  conditions. 

C**********************************************^^^^*****:^*^^***:t^**^^:t^^|^t**** 

SUBROUTINE  rocket ( 

kyOA, yOB.MWA,MWB, ye, yh, temp, pressl,press2,mdot, thrust, Isp.Pjet) 
REAL*8  yOA,yOB,MWA.MWB,ye,yh, 
itemp , press 1 , press2 , mdot , thrust , Isp , P j  et 
REAL’)‘8  MWbar, prat , vex 
DATA  gasR.grav/8. 31439, 9. 8062/ 

MWbar= (yO A*MWA+yOB*MWB ) / ( y h+y e ) 
prat=press2/pressl 

vex=DSQRT(5 . e3*gasR/MWbar*temp*(l-prat*prat*DSQRT(prat) ) ) 

thrust=mdot*vex 

Isp=vex/grav 

P j  et®0 . 5*mdot*vex*vex 

END  ! subroutine  rocket 

C*******************************t******‘**■*************^^********i^■******* 


A. 29  output.for 

c — The  mother  of  all  subroutines.  Basically  has  two  branches:  (1)  if 
c — the  has  selected  howmany=l,  this  subroutine  will  print  anything  of 
c — interest  to  a  file  name  name. out.  where  “name’’  is  the  name  of  the 
c — compound  formed  by  the  two  species  the  user  chose.  The  user  gets  his 
c — original  input  data,  the  electronic,  rotational,  and  vibrational 
c — partition  functions  of  all  the  species,  the  right  hand  sides  of  the 
c — saha  equation,  the  nuclear  fractions  (y’s)  of  the  species,  the  mole 
c — fractions  of  the  species,  the  densities,  the  heats  of  formation  for 
c — the  species,  the  enthalpies,  heat  capacities,  and  entropies,  both  of 
c — the  species  and  for  the  entire  gas,  and  the  rocket  performance.  If 
c — you  W2uit  it  and  its  not  there,  this  program  doesn’t  do  it.  (2)  If  the 
c — user  selects  multiple  temperatures  (howmamy. !=  1),  a  more 
c — convenient,  tabulated  form  of  data  is  outputted  to  several  files, 
c — one  file  per  data  type  per  element.  However,  due  to  the  MS-DOS 
c — limitations  on  the  number  of  possible  files  open  at  one  time,  the 
c — user  must  choose  between  density  and  mole  fraction.  I  chose  these 
c — since  it’s  not  likely  the  user  will  want  both.  Depending  on  the 
c — user’s  choice  of  howmany,  subSUBROUTINEs  open  the  appropriate  files 
c — and  make  the  appropriate  headers . 

(;;tc:<c**:f<*:tc:ti;fltc**4<4:***4<*«*’»:**>«:4t***«4'’t‘:t"t‘******************************4:***** 

SUBROUTINE  output ( 

&yOA , yOB , tempi , temp2 , howmany , pressl , press2 ,mdot , temp , itemp , 

&chA , ZA , rhsA , yA , enthA , CpA , entr A , 

&chB , ZB , rhsB , yB , enthB , CpB , entrB , 

ftname , ZAB , rhsAB , yAB , enthAB , CpAB , entr AB , enthe , Cpe , entre , 
ftenthal .heatcp , entrop ,mf A ,mf B ,mf AB , denA , denB , denAB , 

&dir A , d irB , d ir AB , Pabs , F j  et , thrust , I sp , densmol ) 

CHARACTER  chA*2,chB*2,name*4 
INTEGER  itemp, densmol, howmany 
REAL*8 

&yOA , yOB , tempi , t emp2 , pressl , press2 , mdot , temp , 

&ZA(7,3) ,rhsA(7) ,yA(7) ,enthA(7) ,CpA(7) ,entrA(7) , 

&ZB(7,3) .rhsB(7) ,yB(7) ,enthB(7) ,CpB(7) .entrB(7) , 

&ZAB(2,3) ,rhsAB(2) .yAB(2) ,enthAB(2) ,CpAB(2) .entrAB(2) , 
fcenthe, Cpe, entre, 

4enthal,heatcp, entrop, mfA(7) ,mfB(7) ,mfAB(2) , 

&denA(7) ,denB(7) ,denAB(2) .dirA(7) ,dirB(7) .dirAB(2) , 

&Pabs . P j  et , thrust , Isp 
CHARACTER  filenarae*12 
IMPLICIT  REAL*8  (a-h,l-z) 

DATA  boltzk/1.38066E-23/ 
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INCLUDE  format. for 

c - re-evaluate  the  heavy  and  electron  fractions,  just  to  be  safe. 

CALL  getyeCyA.yeA) 

CALL  getye(yB,yeB) 

CALL  getyh(yA,yhA) 

CALL  getyh(yB,yhB) 
ye=yeA+yeB+yAB(2) 
yh=yhA+yhB+yAB ( 1 ) +yAB (2) 
beta=ye/(yh+ye) 

pktspressl/ (bolt2k*temp* (ye+yh) ) 

c - send  a  teaser  to  the  screen 

WRITE(*, 192) ’  Pabs,  Pjet  (W)  '.Pabs,  Pjet 
WRITE(*,192) '  Thrust(N) , IspCs) ’ .thrust, Isp 
IF  (howmany.EQ.l)  THEN 

c - output  for  single  input 

IF  (itemp.EQ.l)  THEN 

OPEN  (UNIT=90,  FILE^filenaaeCneune, ' .out')) 

END  IF 

c - write  input  info 

WRITE(90 , ' ( A) ' ) ’ Us ing  ' //chA// ' . ' //chB// ' . ' //name 
WRITE(90,*) ’Stoichiometric  Proportions  of  ’//chA//’ , ’//chB 
WRITE(90.’(”  [”,1PE11.3,”.’MPE11.3,”]  ”)’)yOA,yOB 
WRITE(90,*) 'chamber,  exit  pressures  (atm)’ 

WRITE(90,’(”  [”,1PE11.3,”,”,1PE11.3,  ”]  ”) ’)pressl  ,press2 
WRITE(90,*) 'Temperature  range;  templ(K) ,temp2(K) ,howmany ’ 
WRITE(90,'(”  C”,lPE11.3,”,”,lPEil.3,”,",I5,"]  ”)’) 
k  tempi ,temp2.howmany 

WRITE(90,*) 'mass  flow  rate  (kg/s)’ 

WRITE(90,’(”  [”,1PE11.3,”]  ”)’)mdot 

c - pretty  display  of  partition  function 

CALL  Zlshow(chA,ZA) 

CALL  Zlshow(chB.ZB) 

CALL  Z2show(name,ZAB) 

c - show  the  information  of  a  species,  described  by  the  first 

c - parameter  string 

CALL  specshow( ’Right-hand-sides  of  Saha  Equation’, 
k  chA,chB,naune,rhsA,rhsB,rhsAB) 

CALL  specshow( 'Nuclear  fractions, y’ , 
k  chA, chB, name, yA,yB,yAB) 

WRITE(90 , 194) ’ yh , ye , yh+ye , beta ' , yh , ye ,yh+ye , beta 
CALL  specshow( 'Mole  fractions ,y/(yh+ye) ’ , 
k  chA, chB, name, mfA,mfB,mfAB) 
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WRITE(90 , 191) ’ e- ‘ . ye/ (yh+ye) 

CALL  specshowC ’Species  densities  Cd*-3)’, 
ft  chA.chB.name.denA.denB.denAB) 
WRITE(90.191)'e-».pkt*ye 
CALL  specshowC’Heat  of  formation  (J/mol)', 
ft  chA,chB,name.dirA,dirB.dirAB) 

CALL  specshovC' Enthalpies  (J/mol)’, 
ft  chA , chB , name . enthA , enthB , enthAB) 
WRITE(90.191)’e-’.enthe 
WRITE(90,191)name//'  gas’.enthal 
CALL  specshowC’Heat  capacities  (C_p)  (J/K/mol)’, 
ft  chA.chB.name.CpA.CpB.CpAB) 

WRITE(90.191)’e-’,Cpe 
WRITE(90,191)name//’  gas’.heatcp 
CALL  specshowC ’Entropy  diff  (from  STP)  (J/K/mol)’, 
ft  chA , chB , name , Cp A , CpB , CpAB ) 

WRITE(90.191)’e-’.entre 
WRITE(90,191)name//’  gas’.entrop 
WRITE(90,*) 

WRITE(90 . 192) ’Thrust (N) . I8p(s) ’ .thrust , lap 
WRITE(90,192)’Pabs.  Pjet  (W) ’ .Pabs.Pjet 
ELSE 


c - output  for  multiple  inputs 

IF  (itemp.EQ.l)  THEN 

c - open  all  the  appropriate  files  for  output 

CALL  opunhdr ( chA , chB , name , densmol) 

ENDIF 

c - write  all  the  data  to  all  the  files 


WRITE(61 , 139) temp , enthA , enthe 
WRITE(62, 139)temp, enthB, enthe 
WRITE(63, 134) temp, enthAB, enthe 
WRITE(71. 139) temp, CpA.Cpe 
WRITE(72.139)temp.CpB.Cpe 
WRITE(73,134)temp,CpAB.Cpe 
WRITE(81 , 139) temp, entrA.entre 
WRITE(82 , 139) temp , entrB , entre 
WRITE(83 , 134) temp , entrAB , entre 

c - the  MS-DOS-imposed  choice 

IF  (densmol. EQ.l)  THEN 
WRITE(9l , 139)temp,denA,pkt*ye 
WRITE(92, 139)temp,denB,pkt*ye 
WRITE(93 , 134) temp , denAB ,pkt*ye 
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ELSE 

WRITE(91,139)temp,inf  A,ye/(yh+ye) 
WRITE(92.139)temp.mfB,ye/(yh+ye) 

WRITE(93 , 134) temp , mf AB , ye/ (yb+ye) 

EHDIF 

WRITE(94 , 135) temp , Pabs , P j  et . thrust , Isp 
ENDIF 

END  ! subroutine  output 

Qt********************************************************************* 
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A. 30  opunhdr.for 

c — opens  (“opun*’)  all  the  output  files  for  the  houniany!  =  l  case,  and 

c— writes  the  headers  (“hdr"),  such  as  a  line  describing  what  the  file 

c — contains,  with  units,  and  species  column  headings. 

c — file  unit  key 

c — 60*s*enthalpy 

c — 70's=heat  capcity 

c— 80  ’  8®entropy 

c — 90 ’ s=density/molef raction 

c — +l=element  A 

c — +2=element  B 

c— +3=hetero-dia  molecule  AB 

c — output  file  naming  convention  is  as  follows: 

c— (name)_(description)  .out 

C"Where  (name)  is  the  element  or  hetero-dia  molecule  name 
c — and  (description)  is  {eth,cp,etr,den,mol>  for  enthalpy,  etc. 

Qi^*^^i^^^^,^i^^^^*^|^*^^l^H:^^^l:^^***^^*^t^^****^l^m************************************** 

SUBROUTINE  opunhdr( 

AchA , chB , name , densmol) 

CHARACTER  chA*2,chB*2,name*4, 

&filename*12,pref ix*4(3) ,exten*8(4),header*40(4) 

INTEGER  densmol.i, j.ijunit 

pref ix(l)=chA 

prefix (2) =chB 

pref ix(3)=name 

exten ( 1 ) = ’ _eth . out ’ 

exten(2)='_Cp.out ’ 

exten(3)='_etr .out ' 

header(l)=’Enthalpy  (J/mol)’ 

header (2)= ’Heat  capacity,  C_p  (J/K/mol)' 

header (3)=’ Entropy  difference  from  STP  (J/K/mol)’ 

IF  (densmol. EQ.l)  THEN 
exten(4)=’_den.out ' 
header(4)=’Partial  density  (m‘-3)’ 

ELSE 

exten(4)=’_mol . out' 
header(4)=’Mole  fraction’ 

ENDIF 

DO  10  i=l,4 
DO  10  j=l,3 

ijunit=60+10*(i-l)+j 

OPEN  (UNIT=i j unit , FILE=f ilename(pref ix ( j ) , exten(i) ) ) 
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10  CONTINUE 

OPEN  (UNn-94,  FILE=f  ilenameCname,  *_rkt.out')) 

DO  20  i*1.4 
DO  20  j«l,3 

i  j  ujait=60+10*  C  i -1 )  +  j 

20  MRITECijunit ,*)prefix(j),header(i) 

VRITE(94,*)iiaiBe, '  rocket  performance  parameters’ 

991  F0RMAT(2x,’temp(K)',4x.A,9x,A,»+’.8x,A.’+2’,7x.A/+3’.7x.A, 

A’+4’,7x,A. ‘2’ .8x,A. *2+’ ,7x, ’e-') 

993  F0RMAT(2x. »temp(K) » .4x.A.7x,A. ’♦» .6x. ’e-’) 

DO  30  i=l,4 

ijunit=60+10*(i-l) 

WRITECi junit+1 , 991 ) chA , chA , chA , chA . chA , chA . chA 
WRITEC i juait+2 . 991 ) chB , chB , chB , chB . chB . chB . chB 
WRITECi j  unit+3 , 993) name , name 
30  CONTINUE 

994  F0RMAT<2x. ’temp(K) ’ .4x, ’Pabs(W) ’ ,4x, »Pjet(W) ‘ .4x. ’thrust(N) ' . 
ft2x, ’Isp(s) ') 

WRITE(94,994) 

END  'opunhdr 

C***************************************0******t*:^**^i***t****:********** 


A.31  fllename.for 

c — takes  two  strings  (of  arbitrary  length)  and  concatenates  them,  taking 
c— the  spaces  out  of  the  first  string.  This  is  particularly  useful  when 
c—aaking  a  filename  for  the  OPEN  command  that  MS-DOS  will  recognize, 
c — If  there  are  spaces  in  between  the  prefix  and  the  extension,  MS-DOS 
c — makes  a  file  with  no  extension... 


CHARACTER*12  FUNCTION  f ilename(name.exten) 
CHARACTER  name*(*) ,exten*(*) .temp*l2 
INTEGER  ilast 
temp^name 

ilast=INDEX(name, *  ’) 

IF  (ilast. eq.O)  THEN 

temp(LEN(name)+l; )=exten 
ELSE 

temp(ilast : )=exten 
END  IF 

f ilename=temp 
END 
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A.32  specshow.for 

c — displaya  the  contents  of  a  species  variable  to  the  howmany=l  file, 
c— but  in  an  user-interpretable  fashion. 

SUBROUTINE  specshowC 
tstr , chA , chB , name , yA . yB . yAB ) 

CHARACTER  str*(*) ,chA*2,chB*2,nameM 
REAL*8  yA(7),yBC7).yAB(2) 

100  FORMAT  (4x,5(’atomic, ’ ,4x) ,2(*diatomic, ’ ,2x)) 

101  FORMAT  (4x, 'neutral', 4x.4(’  +  M1.9x). ’neutral’. 4x,'+l’) 

102  FORMAT  (A, 1P7(E11 .3)) 

103  FORMAT  (A.53x. 1P2(E11 .3)) 

WRITEOO,*) 

WRITE(90.*)str 

WRITEOO.lOO) 

WRITE(90,101)1,2,3.4 
WRITE(90,102)chA,yA 
WRITE(90,102)chB.yB 
HRITE(90 , 103) name , yAB 
END 

C**************************************************iti**7^*^**^*^**:^*^*i^:^^ 
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A.33  zlshow.for 


c — shows  the  partition  function  for  an  element  in  an  intelligible  way. 

C**************************************************^^*'^*^****^ni*^*0^**m 

SUBROUTINE  Zlshow( 
ftch.Zl) 

CHARACTER  ch*2,partstr*3(3) 

REAL*8  Zl(7.3) 

INTEGER  i 

DATA  partstr/’el  ' . ’vib’ , ’rot’/ 

100  FORMAT  (4x,5(’atomic, ’ ,4x) ,2(’diatomic, * ,2x)) 

101  FORMAT  (4x, ’neutral’ ,4x,4(’+’ ,11, 9x) , ’neutral’ ,4x, '+1 ') 
WRITEOO,*) 

WRITE(90,*) ’Partition  Function:  ’,ch 
WRITE(90,100) 

WRITE(90,101)1,2,3,4 
DO  10  i=l,3 

WRITE(90,’(A,1PE10.3,1P6(E11.3))’) 
ft  partstr(i),Zl(l,i),Zl(2,i),Zl(3,i), 
ft  Zl(4,i),Zl(6,i),Zl(6,i),Zl(7,i) 

10  CONTINUE 
END 
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A. 34  z2show.for 


c — shows  the  hetero-diatomic  species  partition  function;  for  use  with 
c— the  howmany*!  output  file. 

SUBROUTINE  Z2show( 

|[name,Z2) 

CHARACTER  name*4,part8tr*3(3) 

REAL*8  22(7,3) 

INTEGER  i 

DATA  peurtstr/’el  ’ ,  *vib’ , 'rot  V 

100  FORMAT  (4x,5(’atomic, ’ ,4x) ,2(’diatoinic. ' ,2x)) 

101  FORMAT  (4x, 'neutral' ,4x,4( ’+’ ,11, 9x), 'neutral' ,4x, ’+1’ ) 
WRITEOO,*) 

WRITEOO,*) 'Partition  Function:  '.name 
WRITE(90,100) 

WRITE(90, 101)1,2.3,4 
DO  10  i=l,3 

WRITEOO,' (A, 54x.lP2(Ell. 3)) ')part8tr(i),Z2(l,i).  22(2,  i) 

10  CONTINUE 
END 

C********************************************************************** 
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A.35  shutdown.for 


c — shuts  down  the  program,  closing  all  the  open  files,  telling  the  user 
c — which  files  to  look  for,  and  waiting  for  an  <Enter>  before  returning 
c — control  to  DOS.  This  helps  when  running  the  program  from  a  DOS 
c — shell,  like  DCOM,  which  may  automatically  go  back  to  the  interface 
c — screen  and  prevent  the  user  from  seeing  which  files  he  needs  to  see. 
Q*:^m******************************************************************* 

SUBROUTINE  shutdown ( 
fthowmany , chA , chB .name , densmol) 

CHARACTER  chA*2,chB*2,name*4, 
ftfilename*12,pref ix*4(3) ,exten*8(4) ,pause*l 
INTEGER  densmol, i,j,howmany 
IF  (howmany.EQ.l)  THEN 

WRITEC*, ♦) ’Look  for  output  in  file  ( ’ ,f ilenameCname, ’ .out ’ ) , ’ ) ’ 
ELSE 

prefix(l)=chA 
prefix (2) =chB 
prefix(3)=name 
exten(l)=’_eth.out’ 
exten(2)*’_Cp.out ’ 
exten(3)=’_etr .out ’ 

IF  (densmol. EQ.l)  THEN 
exten(4)= ’ _den . out ' 

ELSE 

exten(4)=’_mol . out’ 

ENDIF 

WRITEC*,*) ’Look  for  output  in  files:' 

DO  10  i=l,4 
DO  10  j=l,3 

WRITE(*,*)  ’  ( ’  ,filename(prefix(j)  .extend)) ,  ’)  ’ 

10  CONTINUE 

WRITEC*,*) ’ ( ’ .filename (name, ’_rkt .out’) , ’ ) ’ 

ENDIF 

CL0SE(UNIT=61) 

CL0SE(UNIT=62) 

CL0SE(UNIT=63) 

CL0SE(UNIT=71) 

CL0SE(UNIT=72) 

CL0SE(UNIT=73) 

CL0SE(UNIT=81) 

CL0SE(UNIT=82) 

CLOSE ( UN IT=83) 
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CL0SE(UNIT=90) 

CL0SE(UNIT»91) 

CL0SE(UNIT=92) 

CL0SE(UNIT=93) 

WRITE(*,*) ’Press  <Enter>  to  exit  program...’ 

READ(*,»(A)’)pause 
END  {subroutine  shutdown 

C #  41  *  4i  4:  << ’4' «  *  >tc  *  i)t  *  *  ><1  i|E  Ik  *  *  *  4i  41 4>  *  *  *  * ’tc  «  4c « «  4:  #  *  4:  ♦  * 


100 


• 

A.36 

format. for 

102 

FORMAT 

(214) 

111 

FORMAT 

(A. 14) 

" 

112 

FORMAT 

(A, 214) 

113 

FORMAT 

(A. 314) 

131 

FORMAT 

(1P1(E11.3)) 

132 

FORMAT 

(1P2(E11.3)) 

133 

FORMAT 

C1P3(E11.3)) 

134 

FORMAT 

C1P4(E11.3)) 

135 

FORMAT 

(1P5(E11.3)) 

136 

FORMAT 

(1P6(E11.3)) 

137 

FORMAT 

(1P7(E11.3)) 

138 

FORMAT 

(1P8(E11.3)) 

139 

FORMAT 

(1P9(E11.3)) 

1310 

FORMAT 

(1P10(E11.3)) 

1311 

FORMAT 

(1P11(E11.3)) 

1312 

FORMAT 

(1P12(E11.3)) 

1313 

FORMAT 

(1P13(E11.3)) 

1314 

FORMAT 

(1P14(E11.3)) 

1315 

FORMAT 

(1P15(E11.3)) 

1316 

FORMAT 

(1P16(E11.3)) 

• 

1317 

FORMAT 

(1P17(E11.3)) 

1318 

FORMAT 

(1P17(E11.3)) 

191 

FORMAT 

(A.1PE11.3) 

192 

FORMAT 

(A,1P2(E11.3)) 

193 

FORMAT 

(A,1P3(E11.3)) 

194 

FORMAT 

(A.1P4(E11.3)) 

195 

FORMAT 

(A,1P5(E11.3)) 

196 

FORMAT 

(A,1P6(E11.3)) 

197 

FORMAT 

(A,1P7(E11.3)) 

198 

FORMAT 

(A,1P8(E11.3)) 

199 

FORMAT 

(A.1P9(E11.3)) 

1910 

FORMAT 

(A,1P10(E11.3)) 

1911 

FORMAT 

(A,lPil(Ell.3)) 

1912 

FORMAT 

(A,1P12(E11.3)) 

1913 

FORMAT 

(A,1P13(E11.3)) 

1914 

FORMAT 

(A,1P14(E11.3)) 

1915 

FORMAT 

(A,1P15(E11.3)) 

• 

1916 

FORMAT 

(A,1P16(E11.3)) 

• 

101 

B  Input  Files 


i 


i 


102 


ist  MW(ama)  uD(eV)  ul(eV) 

C  5  12.011  6.21  12.1 

H  2  1.0079  4.47  15.4 

N  4  14.0067  9.75  15.5 

0  5  15.9994  5.11  12.0 

Ar  5  1.800E+01  -5.000E+01  5.000E+01 


Tr,0(K)  Tv,0(K)  Tr,l(K)  Tv,l(K) 

2.618  2668.  2.386  1942. 

87.55  6332.  43.45  3339. 

2.875  3393.  2.779  3175. 

2.079  2273.  2.433  2740. 

l.OOOE+00  l.OOOE+00  l.OOOE+00  l.OOOE+00 


B.2  abxy.in 


uD(eV) 

ul(eV) 

Tr,0(K) 

Tv,0(K) 

Tr.KK) 

Tv,l(K) 

cc 

6.21 

12.1 

2.618 

2668. 

2.386 

1942 

CH 

3.46 

10.6 

20.80 

4112. 

20.39 

3941 

CN 

7.70 

14.1 

2.733 

2976. 

2.728 

2925 

CO 

11.0 

14.0 

2.778 

3121. 

2.844 

3185 

HH 

4.47 

15.4 

87.55 

6332. 

43.45 

3339 

HN 

3.47 

13.6 

24.02 

4722. 

22.08 

4204 

HO 

4.39 

12.9 

27.20 

5377. 

24.16 

4479 

NN 

9.75 

15.5 

2.875 

3393. 

2.779 

3175 

NO 

6.49 

9.26 

2.405 

2739. 

2.873 

3419 

00 

5.11 

12.0 

2.079 

2273. 

2.433 

2740 

ArD 

-5.000E+01 

5.000E+01 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

l.OOOE+00 

CO  b- 


B.3  ar-O.in 


127109.9 
1  0 

5  93143.8 
3  93750.639 
1  94553.707 
3  95399.870 
104102.144 

105462.804 
5  105617.315 

3  106087.305 

S  106237.597 

1  107054.319 

3  107131.755 

5  107289.747 

3  107496.463 

1  108722.668 

1  111667.87 

3  111818.09 

9  112750.22 

7  113020.39 

5  112138.98 

3  114147.75 

5  113426.05 

7  113716.61 

5  114641.04 

7  114821.99 

5  114805.18 

3  115366.90 

5  113468.55 

3  113643.26 

1  114861.67 

3  114975.07 

3  116660.054 

7  116942.815 

5  116999.389 

3  117151.387 

5  117183.654 

1  117563.020 

3  118407.494 

5  118469.117 

3  118459.662 
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1  118870.981 

1  118512.17 

3  118651.447 

9  119023.699 

7  119212.93 

5  118906.665 

3  119847.81 

5  119444.88 

7  119566.11 

5  120619.076 

7  120753.52 

5  120600.944 

3  121011.979 

5  119683.113 

3  119760.22 

I  121096.67 

3  121161.356 

3  120188.34 

S  120188.66 

II  120207.32 

9  120207.77 

7  120229.81 

5  120230.07 

16  120250.15 

16  121653.40 

7  121654.32 

5  121654.58 

3  121068.804 

7  121165.431 

5  121191.92 

3  121257.227 

5  121270.682 

1  121470.304 

3  122609.76 

5  122635.128 

3  122601.290 

1  122790.612 

1  121794.158 

3  121932.908 

9  122036.134 

7  122160.22 

5  122086.974 
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CO  li>  N-  u>  lO  <0 


122514.29 
122282.134 
122329.72 
123505.536 
123557.459 
123372,987 
123815.53 
5  122440.109 

3  122479.459 

1  123873.07 

3  123882.30 

3  122686.20 

5  122686.40 

20  122695.70 

7  122707.94 

5  122708.18 

16  122717.90 

16  124135.74 

7  124137.29 

5  124137.45 

3  123172.09 

7  123205.83 

5  123220.73 

3  123254.99 

5  123261.593 

1  123385.13 

3  124643.54 

5  124658.52 

3  124651.05 

1  124749.89 

1  123508.96 

3  123468.034 

9  123653.238 

7  123773 . 920 

8  123808.60 

5  123826.85 

7  123832.50 

5  125113.48 

7  125150.00 

5  125066.501 

3  125286.28 

5  123903.295 
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3 

123935.97 

1 

125334.75 

3 

125353.31 

3 

124041.20 

5 

124041.38 

20 

124046.64 

7 

124051.44 

5 

124051.65 

16 

124058.36 

16 

125482.70 

7 

125483.16 

5 

125483.34 

3 

124311.72 

7 

124349.04 

5 

124356.73 

3 

124376 . 38 

5 

124381.01 

1 

124439.41 

3 

125783.8 

5 

125791.94 

3 

125777.3 

1 

125831.45 

1 

124526.75 

3 

124554.939 

9 

124609.917 

7 

124649.549 

5 

124603.957 

3 

124788.39 

5 

124692.02 

7 

124715.16 

5 

126064.50 

7 

126089.56 

8 

126053.21 

5 

124771.67 

3 

124782.77 

1 

126202.82 

3 

126211.57 

3 

124857.27 

5 

124857.42 

20 

124860.64 

7 

124865.04 

5 

124865.19 
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lO  CO 


16  124868.77 

16  126294.90 

12  126295.02 

3  125039.60 

7  125054.1 

125059.8 
125072.6 
5  125074.9 

1  125122.54 

4  126524.2 

1  125163.00 

3  125135.898 

9  125219.88 

7  125269.52 

8  125282 . 97 

5  125291.45 

7  125293.65 

5  125329.99 

3  125331 . 93 

8  125386.41 

20  125388.65 

7  125391.04 

5  125391.17 

16  125393.79 

3  125505.5 

12  125519.9 

3  125531.5 

5  125533.8 

1  125561.9 

1  125595.11 

3  125613.12 

-1  -1 
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B.4  ar-l.in 


222820. 

4  0. 

2  1432.0 
2  108722.5 
8  132328.22 
6  132482.12 
4  132631.64 
2  132738.60 
6  134242.62 
4  135086.88 
2  135602.62 
4  138244.51 
2  139259.22 
10  142187.42 
8  142718.01 
6  143108.63 
4  143372.48 
2  144710.90 
4  145669.84 
2  147229.17 
4  147504.12 
6  147876.98 
4  148620.93 
6  148843.29 
8  149180.18 
6  150148.54 
4  150475.82 
6  151088.18 
6  155044.07 
4  155352.04 
2  155709.02 
8  157234.93 
6  157674.30 
4  158168.71 
2  158429.05 
6  158731.2 
4  159394.32 
2  159707.46 
4  160240.35 
4  161049.65 
2  161090.31 


no 


2  167308.66 
6  170401.88 
8  170531.29 
4  172214.74 
2  172817.12 
6  172336.47 
4  172830.63 
4  173348.78 
6  173394.33 
4  174410.74 
2  174821.94 
4  179593.09 
2  179932.83 
6  181595.04 
4  182223.06 
2  182952.14 
4  183091.83 
2  183915.58 
8  183676.42 
6  183798.22 
4  183986.83 
2  184193.12 
2  184094.10 
10  185093.92 
8  185625.47 
6  186074.06 
4  186341.39 
2  186172.32 
4  186471.32 
6  186891.92 
6  186728.28 
4  186750.78 
8  186817.12 
6  187589.62 
2  189935.62 
4  190593.62 
4  190106.84 
2  190196.80 
10  190508.00 
2  191708.46 
4  191975.16 
2  192334.09 


6  192557.77 
4  192712.93 
10  194800.97 
8  194822.95 
6  194862.31 
4  194997.65 
8  194883.96 
6  195032.13 
4  195298.62 
2  195282.50 
6  195865.61 
4  195867.73 
4  196077.40 
2  196091.04 
4  196622.78 
6  196633.93 
8  198595.91 
10  198604.78 
6  198813.17 
4  199138.92 
2  200111.16 
2  199447.56 
4  199982.96 
4  199525.96 
6  199680.58 
4  200032.65 
2  200624.00 
8  200139.84 
6  200235.70 
4  204418.50 
2  204515.81 
10  204586.40 
2  205243.96 
6  208592.90 
4  212932.88 
6  212934.30 
-1  -1 
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329965.80 

F  0. 

3  1112.40 
1  1570.20 
5  14010. 

1  33267. 

5  113800.70 
3  114797.60 
1  115328.40 

3  144023. 

4  144882.93 

5  144885.97 
7  144892.95 
9  144907.00 
7  156917.62 
5  156924.68 
3  157031.40 
5  174375.00 
3  180679.00 
9  186402.15 
7  186657.20 
5  186903.05 
3  187171.12 
5  187823.05 
7  188714.05 
9  188517.32 
3  196589.20 
5  196613.91 
7  196679.80 
3  204563.53 
5  204649.24 
7  204797.37 
3  204727.47 
5  207233.09 
3  207532.15 
1  207673.16 
5  209151.82 
3  209127.04 
1  209166.35 
7  210212.26 
5  211004.85 


113 


3  211563.83 
5  213950.87 
3  214346.70 
1  214568.49 
3  225155.18 
5  225147.93 
7  225402.59 
5  226355.96 
7  226503.22 
9  226646.06 
5  231341.80 
3  231627.30 
1  231754.80 
3  239193.48 
3  240150.66 
5  240257.59 
7  240291.66 

I  242923.96 

3  243145.76 
5  243424.97 

4  246029.76 

5  246033.79 
7  246036.64 
9  246046.57 
5  250712.27 
3  252272.92 
5  252253.69 
7  252289.02 
3  252575.88 
5  266722.80 
7  266877.50 
9  267071.22 
7  267782.10 
9  267833.20 

II  267895.82 
3  268978.80 
5  269012.80 
7  269000.80 
5  271507.88 
3  271672.08 
1  271696.22 
3  272068.45 
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3  272127.82 
5  272188.16 
7  272250.90 
5  281461.97 
16  281473.82 
1  281947.88 
3  282000.26 
5  282099.14 
7  283919.78 
5  284096.26 
3  284118.51 
1  285831.20 
3  285882.00 
5  286009.21 
-1  -1 
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B.6  ar-3.in 


482400. 

4  0. 

4  21090. 

6  21219. 

2  34854. 

4  35035. 

6  117564. 

4  118515. 

2  119044. 

4  145921. 

6  146000. 

4  166356. 

2  167444. 

2  177833. 

2  250219.45 
4  250906.60 
6  251972.00 
2  256093.29 
4  257348.89 
6  268151.38 
4  268171.38 
2  285960.17 
4  286228.80 
6  286751.68 
8  287555.83 
2  289125.88 
4  289237.82 
6  289834.68 
4  290256.45 
6  291667.73 
4  291748.70 
2  295674.54 
4  295806.77 
2  299563.20 
6  304074.29 
8  304399.90 
6  306236.28 
4  306308.25 
-1  -1 
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B.7  ai:-4.in 


605100. 

1  0. 

3  765. 

5  2032. 

5  16301. 

3  121632. 
5  121678. 
7  121840. 
5  141764. 

4  141773. 
3  191537. 
3  195356. 

5  217578. 
3  218286. 
1  218642, 
3  224216. 
5  224505. 
7  224717. 
1  295742. 
3  296249. 
5  297893. 
3  301300. 
-1  -1 


117 


B.8  c-O.in 


90878.3 
1  0.0 
3  16.4 
5  43. S 
5  10193.70 
1  21648.4 
5  33735.2 
1  60333.80 
3  60353.00 
5  60393.52 
3  61982.20 
7  64088.56 
5  64093.19 
3  64092.01 
3  68858. 

3  69689.79 
5  69710.99 
7  69744.40 
3  70744.26 
1  71352.81 
3  71365.23 
5  71385.70 
5  72611.06 
1  73976.23 
9  75256.3 
5  77680.5 
1  78105.23 
3  78117.06 
5  78148.36 
5  78199.34 
7  78215.82 
9  78250.22 
3  78300.8 
5  78307. 

7  78316, 

3  78338. 

7  78531. 

3  78727.91 
5  79311.10 
3  79319.06 
1  79323,32 


118 


3  80173.29 
5  80192.49 
7  80222.74 
3  80563.57 
3  81105.70 
1  81311.52 
3  81326.33 
5  81344.48 
S  81770.36 
1  82252.31 
5  83500. 

7  83761. 

3  83830. 

5  83837. 

7  83847. 

3  83882.5 
7  83949. 

3  84032. 

5  84102.6 
3  84112. 

3  84852.13 
5  84952. 

7  84986.2 
5  85400.38 
1  85625.84 
5  86187. 

5  86319. 

7  86326.7 
5  86371.3 
7  86396. 

3  86413.96 
7  86450. 

3  86491. 

5  86504, 

3  86517. 

5  87632. 

5  87706. 

7  87713. 

5  87752. 

7  87773. 

3  87795.3 
7  87807. 


5  87830. 

3  87839. 

3  87831.3 
5  88541.8 
7  88547. 

7  88607. 

7  88624. 

3  88632.44 
5  88639. 

7  89081. 

5  89082. 

7  89146. 

7  89155. 

5  89158. 

5  89450. 

7  89514. 

7  89517. 

7  89779. 

7  89968.4 
-1  -1 


i 
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B.9  c-l.in 

196659.0 
2  0.0 
4  64.0 

1  46000.2 
4  43021.8 
6  43050.7 
6  74930.9 
4  74933.2 

2  96494.1 
2  110625.1 
4  110666.3 
2  116537.88 
2  131724.68 
4  131735.81 
4  142024.4 
4  145549.99 
6  145551.44 
6  150462.8 
4  150467.9 
2  157234.43 
2  162518.70 
4  162524.62 
2  166964.70 
4  166988.46 
6  167033.43 
4  168123.92 
6  168124.33 
2  168731.6 
4  168750.2 
14  168979.05 
2  173348.18 
2  175287.9 

4  175295.2 
2  178194.1 
4  178220.8 
6  178494.8 
14  178956.46 
2  181258. 

2  181694.50 
4  181709.20 
6  181734.21 


121 


8  181770.48 
2  182025.0 
4  182044.5 
6  184064.9 
14  184376.20 
4  184688.69 
2  186425.02 
4  186441.32 
6  186463.75 
4  188579.3 
6  188612.7 
2  194571.9 
4  195750.8 
6  195765.1 
8  195784.7 
10  195812.3 
2  196556.2 
4  196561.8 
6  196570.5 
8  196580.8 
-1  -1 
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B.IO  c-2.iri 


386159.7 
1  0.0 
1  52315.0 
3  52338.0 
5  52394.8 
3  102351.4 
1  137374.0 
3  137403.4 
5  137450.5 
5  145875.1 
1  182520.2 
3  238160.7 
1  247169.5 
3  258931.4 
1  259653.8 
3  259659.3 
5  259672.1 
3  269957.6 
5  269959.7 
7  269962.9 
5  276482.7 
1  308162.9 
3  308196.2 
5  308264.8 
3  309404.5 

3  310005.2 
1  311720.7 

4  317743. 

5  317748. 

3  319719.4 
3  321358.8 
5  321375.1 
7  321398.6 
5  321949.1 
7  321955.8 
9  321964.7 
3  322403.1 
7  322701.1 
3  323024.0 
5  323049.4 
7  323088.2 


123 


5  324212.0 
3  327225.7 
1  329633.1 
3  329654.2 
5  329690.9 
5  332690.3 
5  333116.4 
5  333333.4 
7  333358.4 
9  333395.0 
3  337602.9 
5  337616.4 
7  337636.7 
3  339881. 

5  340049.5 
3  340075.8 
1  340090.3 
7  341368.5 
3  343255.7 
5  344181. 

I  345093.9 
7  345444. 

9  346525.1 

II  346253.0 
9  346577.5 
5  346656.0 
3  346713.1 
5  347099.5 
7  347101.3 
9  347103.7 
7  348859.5 
3  354796. 

3  357088. 

7  358046. 

9  358638.3 
11  358639.0 
9  358688.9 
5  358725.5 
9  358800. 

7  359122.2 
3  363561. 

3  364896. 


124 


15  365585. 
5  366027.0 
3  369926. 
15  370438. 
15  373748. 
5  376637. 

3  381104.8 
5  381919. 

7  381958. 

3  384313. 

5  384350. 

5  385637.5 
5  385816.2 
-1  -1 


B.ll  c>3.in 


520177.8 
2  0.0 
2  64484.2 
4  64591.3 
2  302847.9 
2  320048.5 
4  320080.0 
4  324880.2 
6  324890.9 
2  401346.7 
2  408308.9 
4  408322.2 
4  410333.8 
6  410338.2 
8  410434.1 
2  445366.1 
2  448854. 

4  448861. 

6  449887.4 
14  449938.2 
18  449948.4 
2  468765. 

6  470763. 

10  471368. 
14  471403.0 
18  471407.4 
22  471407.9 
2  482659. 

6  483931. 

10  484309. 
14  484343.8 
18  484346.6 
22  484346.9 
6  492473. 

14  492743. 
36  492745. 
-1  -1 


126 


B.12  c-4.in 


3162450. 

1  0. 

3  2411266. 

1  2455165. 

3  2455152. 

5  2455288. 

3  2483240. 
15  2857308. 
3  2859350. 

3  2991680. 

3  3053060. 

3  3086420. 

3  3106750. 

3  3118760. 
-1  -1 


127 


B.13  h-O.in 

109678.758 
2  0.000 
8  82259 
14  97492.3 
22  102823.9 
32  105291.6 
44  106632.15 
58  107440.42 
74  107965.03 
-1  -1 


128 


B.14  h>l.in 

100000000.000 
1  0.000 
-1  -1 


129 


B.15  n-O.in 

117345 
4  0.000 

6  19223 
4  19231 
6  28840 
2  83285.5 
4  83319.3 
6  83366 
2  86131.4 
4  86223.2 
6  88109.5 
4  88153.4 
2  88173 
2  93582.3 
2  94772.2 
4  94794.8 
6  94832.1 
8  94883.1 
2  95476.5 
4  95494.9 
6  95533.2 
4  96751.7 
4  96788.2 
6  96864.2 
2  97770.1 
4  97805.8 
6  99665 
4  99658 
2  103618.1 
4  103668.1 
6  103736.8 
2  104142.2 
4  104227.4 
4  104615.4 
2  104654.9 
4  104665 
6  104684 
8  104718 
10  104767 
6  104810.9 
8  104882.7 


130 


2  104864 
4  104890 
6  104957 
2  104987 
4  104998 
6  105011 
8  105020 
4  105120.8 
6  105144.3 
2  106478.6 
2  106760.5 
4  106780.1 
6  106816.1 
8  106870.7 
2  106982,7 
4  106998.3 
6  107039 
4  107447.2 
2  109813.5 
4  109857.8 
6  109927.9 
2  110029.2 
4  110108.5 
4  110196 
6  110214 
8  110248 
10  110304 
2  110221 
4  110275 
6  110288 
8  110339 
4  110221.7 
2  110244.6 
6  110311 
8  110373 
2  110325 
4  110351 
6  110403 
4  110448.3 
6  110470.5 
4  110521.9 
6  110545.8 


131 


2  112294.8 
4  112320.8 
2  11256S.9 
4  112610.6 
6  112682.6 
2  112735 
4  112823 
4  112751 
6  112763 
8  112799 
10  112862 
4  112801 
2  112816 
6  112820 
8  112890.2 
6  112825 
6  112825 
8  112892 
2  112855 
4  112874 
6  112912 
4  112929.2 
6  112947.5 
2  114015 
4  114072 
6  114146 
2  114130 
4  114163 
28  114160 
12  114182 
8  114248 
4  114193 
2  114209 
6  114196 
8  114275 
4  114232.2 
6  114290.5 
6  114259 
6  114274 
2  114809 
4  114890 
6  114942 


132 


6  114950 
20  114988 
14  115004 
-1  -1 


B.16  n-l.in 


238846.7 
1  0.000 
3  49.1 
5  131.3 
5  16315.7 
1  32687.1 
5  47167.7 

7  92237.9 
5  92251,3 
3  92252.9 

8  109218.2 
1  109224.8 
5  144189.1 
1  148909.37 
3  148940.97 
5  149077,33 
3  149188.74 
3  155129.9 
3  164611.6 
3  166522.48 
S  166583.26 
7  166679.45 
3  166765.7 
3  168893.04 
1  170573.38 
3  170608.63 
5  170667 

5  174212.93 
1  178274.17 
5  186512.38 
7  186571.8 

9  186653.35 
5  187092.2 
3  187438.34 
5  187462.38 
7  187492.72 
5  188858.09 
3  188909.89 
1  188937.95 
7  189336 

3  190121.15 


134 


1  196541.09 
3  196592.88 
5  196712.17 
3  197859.28 
3  202169.9 
3  202714.94 
5  202765.86 
7  202862.06 
1  203164.7 
3  203188.8 
5  203259.7 
3  203532.8 
5  205350.7 
3  205982.1 
5  206038.1 
7  206108.7 
1  206327.5 
5  209675.3 
7  209739.5 
9  209825.3 
5  209926.92 
3  210239.8 
5  210266.3 
7  210301.9 
5  210705.4 
3  210T51 .5 

I  210777 

7  211030.9 
5  211033.71 
7  211057.07 
9  211061.03 
7  211104.8 
7  211288.02 
9  211295.65 

II  211390.77 
3  211335.5 

9  211402.89 
7  211411.25 
5  211416.2 
3  211487.28 
5  211491.16 
1  211750.2 


3  211780.6 
5  211828.8 

I  214212.4 
3  214258.2 
5  214385.3 
3  214828 
15  220717 
12  221070.2 
9  221074.3 
7  221137.6 
7  221227.7 
9  221232.7 

II  221302.2 
9  221312.1 

1  224027.1 
3  224042.9 
5  224072.3 
7  224115.4 
9  224169.3 
3  225987.1 
5  226011.2 
7  226055.2 
5  230223 
-1  -1 


136 


B.17  n>2.in 


382625.5 
2  0.000 
4  174.5 
2  57192.1 
4  57252 
6  57333.2 
6  101023.8 
4  101031.5 
2  131003.5 
2  145876.1 
4  145986.5 
4  186802.3 
6  203072.3 
4  203088.9 
2  221302.4 
2  230404.5 
4  230408.6 
2  245665.7 
4  245701.7 
4  267238.5 
6  267244.4 
2  287535.6 
4  287598.1 
6  287713.9 
2  297150.2 
4  297263.1 
2  301088.2 
2  309132.6 
4  309185.8 
2  309662.8 
4  309698.3 
6  309760.5 
8  309856.7 
2  311691.3 
4  311716.1 
4  314224 
2  317299.9 
4  317343.4 
6  317402.3 
4  317750.8 
6  317781.8 


137 


14  320287.5 
4  320977.4 
6  321065.8 
2  327056.8 
4  330238.4 
6  330273.5 
8  330325.3 
10  330396.7 
2  332796.6 
4  332810 
6  332832 
8  332860.3 
2  333713.1 
4  334542.2 
6  334568.9 
6  336213.4 
4  336268 
2  336303.1 
6  339744.4 
8  339855.7 
4  341946.2 
6  341947.9 
4  342693 
2  342763.7 
14  342752 
18  343116 
10  354517 
14  354955.7 
18  355214 
2  368525.6 
4  368588.3 
6  368704.8 
4  373342 
6  373376 
2  374747.4 
4  374805,3 
2  376756.6 
4  376803.3 
6  376863.8 
8  376953.3 
2  377591 
4  377608 


138 


B.18  n-3.in 

624851 
1  0.000 
1  67136.4 
3  67199.6 
5  67343.8 
3  130695 
1  175463.5 
3  175536.7 
5  175661.5 
5  188885 
1  235370 
3  377206 
1  388858 
3  404521 
1  405893.2 
3  405909 
5  405944 
3  419967.8 
5  419971.3 
7  419979.4 
5  429158 
1  465223 
3  465300.6 
5  465463.4 
3  473032 
3  480880 
5  484394 
7  484525 

3  487542 

4  494240 

5  494338 
5  498315 
5  499708 
21  499851 
9  503625 
3  505487 
5  505518 
7  505561 
7  506292 
3  507022 
15  511384 


140 


5  511440 

4  511493 

5  514638 
5  516631 
7  516639 
9  516650 
3  519414 

7  521868 
3  550218 
15  552731 
27  554419 
15  574940 
5  591043 

8  593665 
7  593704 
-1  -1 


B.19  o>0.in 

109836.7 
5  0. 

3  158.5 

1  226.5 

5  15867.7 

1  33792.4 

5  73767.81 

3  76794.69 

3  86625.35 

5  86627.37 

7  86631.07 

5  88630.84 

3  88630.30 

1  88631.00 

5  95476.43 

3  96225.5 

9  97420.24 

7  97420.37  multiple? 

5  97420.50  mult 
7  97488.14  mult 
3  99092.64 
5  99093.31 
7  99094.52 
5  99680.4 
7  101135.04 
5  101147.21 
3  101155.10 
5  102116.21 
3  102411.65 
5  102661.63 
9  102865.09  multiple? 
7  102908.14  multiple? 
5  103869.4  mult 
5  105019.0 
3  105164.90 
9  105385.3  multiple? 

7  105408.58  mult 
5  105911.3 
5  106545.1 
3  106627.9 
9  106751.2 


142 


7  10676S.8 
5  107445.4 
3  107497.1 
9  107573.1 
7  107582.7 
5  108021.4 
3  108057.6 
9  108105.7 
7  108116.6 
5  108412.0 
3  108436.1 
9  108470.2 
7  108477.8 
5  108688.4 
3  108707.3 
9  108731.5 
7  108734.4 
-1  -1 


143 


283550.9 
4  0. 

6  26808.4 
4  26829.4 
4  40466.9 
2  40468.4 
6  119837.7 
4  120001.1 
2  120083.5 
6  165987.7 
4  165996.0 
2  185235.36 
4  185340.68 
6  185499.20 
2  188688.36 
4  189068.37 
2  195710.4 
2  203942.21 
2  206730.80 
4  206786.34 
6  206877.90 
8  207002.52 
6  206971.3 
4  206972.3 
2  208346.17 
4  208392.27 
6  208484.24 
4  211521.98 
6  211712.66 
4  212161.94 
4  212593.2 
2  212762.4 
2  214169.74 
4  214229.48 
2  226851. 

6  228723.3 
8  228746.9 
6  229946.6 
4  229968.2 
4  231296.05 
6  231350.08 


144 


8  231427.99 
10  231530.26 
6  232462.83 
4  232536.06 
2  232602.57 
2  232480.1 
4  232526.7 
2  232711.70 
4  232745.98 
6  232747.51 
8  232753.86 
6  232796.27 
8  232959.26 
4  233430.10 
2  233544.09 
4  234402.48 
6  234454.45 
2  238626.32 
4  238731.54 
6  238892.96 
2  240328.75 
4  240516.28 
6  245395.5 
2  245767.80 
4  245816.29 
6  245902.85 
8  246028.95 
4  248009.1 
6  248185.3 
2  248425.35 
4  248514.23 
6  250251. 

8  251220.9 
6  251224.1 
10  252607.7 
8  252608.9 
4  253046.23 
6  253048.35 
2  253789.51 
4  253791.87 
8  254481.5  multiple? 
10  254590.7 


10  2S4895.2  multiple? 
4  254982.2 
6  2SS104.6 
4  255140.9 
2  255162.6 
4  255172.5 
2  255281.4 
6  255301.3 
8  255465.2 
2  255622.4 
6  255689.6 
4  255812.2 
8  255691.4 
6  255213.1 
4  255913. 

2  255912.0 
6  255755.8 
8  255759.4 
10  255827.6 
12  255977.5 
8  255829.4 
10  255983.6 
4  255843.1 
6  255897.2 
4  256083.5 
6  256087.6 
8  256123.1 
10  256136.2 
6  256125.8 
8  256143.3 
2  257693.7 
4  257797.9 
6  257963.8 
2  258408.6 
4  258601.7 
6  259286.2 
4  259287.0 
4  260959. 

6  261042. 

8  261180. 

4  261261.7 
6  261354.3 


146 


4  26ie9t;S 
6  261869.4 
10  265220.3 
6  265431.5 
6  265468.2 
8  265578. 

8  265639. 

6  265705. 

4  265762. 

2  265859. 

6  265665. 

8  265691. 

10  265761. 
12  265925. 

8  265763.0 
10  265930.2 
6  265856. 

4  265928. 

6  265961. 

8  265985. 

10  265999. 

6  265988. 

8  265999. 

4  267763.39 
6  267770.85 
8  267783.40 
6  274739.2 
8  274782.4 
10  274920. 

6  275611. 

18  275841.3 
14  275879.6 
10  275951. 

2  275997. 

10  276066.3 
22  276109.1 
6  276263.9 
10  278140 
-1  -1 


147 


B.21  0-2. in 


443193.5 

1  0. 

3  113.4 
5  306.8 
5  20271.0 
1  43183.5 
5  60312.1 
7  120025.4 
5  120052.6 
3  120056.5 
5  142381.7 
3  142382.8 
1  142396.9 
5  187049.4 
3  197086.7 
3  210458.5 
1  267257.29 
3  267375.65 
5  267632.59 
3  273080.07 
5  283758.9 
3  283976.6 
1  284073.3 
3  290956.62 
3  293865.26 
5  294001.60 
7  294221.65 
3  297557.50 
5  298289.4 
1  300228.21 
3  300310.31 
5  300440.85 
5  306584.8 
1  313801.07 
5  324462.46 
7  324658.25 
9  324836.41 
5  324734.22 
3  327227.94 
5  327277.18 
7  327350.90 
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5  329467.98 
3  329581.98 
1  329643.43 
7  331820.2 
3  332777.1 
3  338565.87 
5  338690.34 
7  338851.50 
1  343302.6 
1  350026.1 
3  350122.9 
5  350302.3 
1  356732. 

3  356838. 

5  357111. 

3  358667.4 
3  363266.8 
1  365515.76 
3  365550.60 
5  365619.12 
7  365719.16 
9  365846.46 
3  365723.9 
3  366486.91 
5  366594.01 
7  366801.04 
3  367952.20 
3  368526.37 
5  368583.63 
7  368684.75 
1  370326.7 
3  370415.7 
5  370524.2 
5  370900.6 
1  373046.2 
3  374575. 

5  374662.5 
7  374798.6 
5  376067.66 
5  377375. 

5  377687. 

5  378408.5 


149 


3  378420.9 

I  378438.1 
3  379232. 

5  379293. 

7  379356. 

5  380706. 

7  380782. 

3  381086. 

5  392221. 

3  392778. 

3  394090. 

5  394126. 

7  394195. 

3  394516.45 
5  394555.15 
7  394612.70 
9  394688.44 

II  394780.47 
1  398135.0 

3  398131.4 
5  398127.3 
7  398137.4 
9  398218.8 
7  398474.3 
5  398544.3 
3  398582.8 
5  400354.8 
3  400464.7 
1  400518.4 
5  401379. 

7  401475.4 
9  401609.1 
5  401530. 

5  401787. 

7  402530. 

7  403374. 

3  403526. 

3  405805.1 
5  405834.1 
7  405833.0 
5  414675. 

7  415181. 
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7  422977. 

7  424998. 

5  426338. 

3  428487. 

5  428606. 

7  428769. 

3  430025. 

3  437015.0 
3  438241.0 
5  438303.2 
7  438395.2 
9  438517.5 
3  439278.1 
5  439329.5 
7  439427.6 
7  442710. 
-1  -1 


B.22  o>3.in 


624396.5 
2  0.0 
4  386.5 
2  71177.0 
4  71308.4 
6  71492.9 
6  126936.3 
4  126950.3 
2  164366.9 
2  180481.3 
4  180724.6 
4  231275.1 
6  255156.7 
4  255186.0 
2  289016.1 
4  289024.0 
2  357614.8 
2  390161.1 
4  390248.2 
4  419533.5 
6  419550.2 
2  438588.5 
4  438723.6 
6  438970.5 
2  452808.0 
4  453073.0 
2  467231.1 
4  467346.5 
2  468075.4 
4  468154.2 
6  468289.7 
8  468499,4 
4  474217.8 
2  478587.7 
4  478682.2 
6  478811.3 
4  482667.5 
6  482923.1 
2  485823.1 
2  492880. 

4  494907.5 


152 


6  494986.3 
8  495098.7 
10  495252.8 
2  499506.4 
4  499535.3 
6  499582.0 
8  499646.6 
4  501511.3 
6  501566.4 
6  503834.5 
4  503947.9 
2  504021.7 
4  510560. 

6  510567. 

6  510746.1 
8  510978.5 
4  514217. 

2  514368. 

2  518684. 

4  518690. 

2  539368. 

4  547311. 

6  547336. 

2  549792. 

4  549855. 

6  552034. 
14  552490 
2  554461 
2  568638. 

4  568773. 

6  569020. 

8  570791, 

2  573696. 

4  573907. 

6  574373. 

2  575204. 

4  575373. 

4  575819. 

6  575853. 

2  576591. 

4  576735. 

6  576947. 


153 


2  581721. 

4  581743. 

4  584552. 

6  584768. 
14  587850. 
2  590071. 

8  591767. 

6  592999. 

4  593627. 

6  593708. 

6  594007. 

8  594080. 

4  594337. 

6  594542. 

6  596299. 

8  596477. 

2  597254. 
14  597352. 
4  597726. 

2  597863. 

4  600092. 

6  600106. 

8  602977. 

6  602434. 

6  615431. 

4  615460. 

4  616588. 
-1  -1 
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B.23  o-4.itt 


918702. 

1  0. 

1  82121.2 
3  82257.9 
5  82564.1 
3  158798. 

1  213641.7 
3  213797.4 
5  214066.2 
5  231722. 

1  287909, 

3  547150.0 
1  561278. 

3  580826. 

1  582983.6 
3  583019.9 
5  583097.2 
3  600925.5 
5  600936.6 
7  600956.1 
S  612617. 

1  653099.7 
3  653262.2 
5  653605.0 
3  664486. 

3  672695. 

3  677333. 

5  677532. 

7  677847. 

3  684124. 

1  689585.6 
3  689699.6 
5  689890.3 
5  694646. 

5  697170. 

3  704360. 

5  704424. 

7  704527. 

1  707630. 

5  708154. 

3  708296. 
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1  708379. 
7  712967 
3  709277. 
3  722666. 
1  731667. 
3  736108. 
5  736126. 
3  737883. 
3  742401. 
5  742407. 
7  742421. 
5  746280. 
7  749857. 
3  796263. 
3  802452. 
7  806625. 
5  808351. 
3  824280. 
3  829588. 
3  831047. 
5  831213. 
7  831504. 
3  832251. 
3  835151. 
5  835321. 
5  837834. 
5  837864. 
3  839616. 
7  840832. 
7  841220. 
3  841280. 
5  841374. 
7  841497. 
5  842105. 
5  843290. 
3  843397. 
1  843449. 
7  847129. 
3  847465. 
3  860874. 
7  861975. 
5  862419. 
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A 


A 


3  874447. 
7  875365. 
3  898580. 
7  899671. 
5  901344. 
5  902442. 
5  902592. 
7  904497. 
7  906404. 
-1  -1 
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